{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "BBF0E215DDC04ED79CA38FB7C075CD4C",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "# <center> Lecture14 : Hierarchical Models </center>  \n",
        " \n",
        "## <center> Instructor: Dr. Hu Chuan-Peng </center> "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "3788A3AC09C74A8BA90F0D6018724BAE",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "在本节课中，我们将介绍层级模型(hierarchical model)。主要内容：  \n",
        "- 了解层级数据结构的形式  \n",
        "- 三种对待处理层级数据的思路：完全池化(complete pooling),  非池化(no pooling), 部分池化(patial pooling)  \n",
        "- 了解组间变异(between variability)和组内变异(within variability)的差异  \n",
        "- 通过pymc实现三种不同的模型，并理解分层模型的意义  \n",
        "\n",
        "> 注：层级模型有许多别名，多层线性模型（Hierarchical Linear Model，HLM），也叫多水平模型（Multilevel Model，MLM），线性混合模型（Linear Mixed Model）混合效应模型（Mixed Effects Model）随机效应模型（Random Effects Model)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "9EC04AC26DD1410FBC563957AA9BCD10",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## 层级数据结构"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "9DC2E91D65084813B5F945B76150615E",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "在心理学实验中，层级数据或分组数据(hierarchical or grouped data)十分常见，例如：  \n",
        "\n",
        "1. 设计不同的实验条件，在同一实验条件下收集多个数据(实验条件为组)  \n",
        "2. 划分不同的人群，在同类人群中收集多个数据(人群为组)  \n",
        "3. 在一种实验条件内，被试需要重复对同种刺激做出反应(被试为组)  \n",
        "\n",
        "![Image Name](https://cdn.kesci.com/upload/image/rmcuwtjhi2.png?imageView2/0/w/960/h/960)  "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "25D7894F460241DF9D186C7D6E9331EA",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "在我们之前使用的数据中，层级数据表现为来自多个站点的自我控制分数  \n",
        "\n",
        "> * 数据来源: Hu, C.-P. et al. (2018). Raw data from the Human Penguin Project. Open Science Framework. https://doi.org/10.17605/OSF.IO/H52D3   \n",
        "> * 自我控制量表来源：Tangney, J. P., Baumeister, R. F. & Boone, A. L. High self-control predicts good adjustment, less pathology, better grades, and interpersonal success. J. Pers. 72, 271–324 (2004)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "collapsed": false,
        "id": "3834C39D1B804BD98F99AD492E81B23E",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "skip"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n"
          ]
        }
      ],
      "source": [
        "# 导入 pymc 模型包，和 arviz 等分析工具 \n",
        "import pymc as pm\n",
        "import arviz as az\n",
        "import seaborn as sns\n",
        "import scipy.stats as st\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import xarray as xr\n",
        "import pandas as pd\n",
        "import ipywidgets\n",
        "import bambi as bmb\n",
        "\n",
        "# 忽略不必要的警告\n",
        "import warnings\n",
        "warnings.filterwarnings(\"ignore\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "collapsed": false,
        "id": "8C7E5FFDD44B45B3B518173548771D22",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>Site</th>\n",
              "      <th>scontrol</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>Bamberg</td>\n",
              "      <td>41</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>Bamberg</td>\n",
              "      <td>36</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>Bamberg</td>\n",
              "      <td>31</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>Bamberg</td>\n",
              "      <td>40</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>Bamberg</td>\n",
              "      <td>43</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>...</th>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1518</th>\n",
              "      <td>Zurich</td>\n",
              "      <td>45</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1519</th>\n",
              "      <td>Zurich</td>\n",
              "      <td>39</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1520</th>\n",
              "      <td>Zurich</td>\n",
              "      <td>42</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1521</th>\n",
              "      <td>Zurich</td>\n",
              "      <td>42</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1522</th>\n",
              "      <td>Zurich</td>\n",
              "      <td>42</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "<p>1523 rows × 2 columns</p>\n",
              "</div>"
            ],
            "text/plain": [
              "         Site  scontrol\n",
              "0     Bamberg        41\n",
              "1     Bamberg        36\n",
              "2     Bamberg        31\n",
              "3     Bamberg        40\n",
              "4     Bamberg        43\n",
              "...       ...       ...\n",
              "1518   Zurich        45\n",
              "1519   Zurich        39\n",
              "1520   Zurich        42\n",
              "1521   Zurich        42\n",
              "1522   Zurich        42\n",
              "\n",
              "[1523 rows x 2 columns]"
            ]
          },
          "execution_count": 2,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# 通过 pd.read_csv 加载数据 Data_Sum_HPP_Multi_Site_Share.csv\n",
        "try:\n",
        "  df_raw = pd.read_csv('/home/mw/input/bayes20238001/Data_Sum_HPP_Multi_Site_Share.csv')\n",
        "except:\n",
        "  df_raw = pd.read_csv('data/Data_Sum_HPP_Multi_Site_Share.csv')\n",
        "\n",
        "df_raw[[\"Site\",\"scontrol\"]]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "68CE6D58228B499DA91F9699959FCF52",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 数据可视化  \n",
        "\n",
        "将所有站点的自我控制分数进行可视化，可以发现：  \n",
        "\n",
        "* 有的组存在较多极端值  \n",
        "\n",
        "* 组与组之间的均值不同"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "collapsed": false,
        "id": "1A3CB9EB5742499BAF37B435A69310AA",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/77CBDDCB50124DAAB15A59A575A62555/s5jqbxv522.png\">"
            ],
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "sns.boxplot(data=df_raw,\n",
        "            x=\"Site\",\n",
        "            y=\"scontrol\")\n",
        "\n",
        "plt.xticks(rotation=90) \n",
        "sns.despine()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "B4DF0FDFFFBB493291CF3F06BA9B27EC",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "为了方便之后演示，我们仅使用其中5个站点的数据："
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "collapsed": false,
        "id": "A106D4AF4CA24CD2931905A4B0571C35",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th>age</th>\n",
              "      <th>anxiety</th>\n",
              "      <th>anxiety_r</th>\n",
              "      <th>artgluctot</th>\n",
              "      <th>attachhome</th>\n",
              "      <th>attachphone</th>\n",
              "      <th>AvgHumidity</th>\n",
              "      <th>avgtemp</th>\n",
              "      <th>avoidance</th>\n",
              "      <th>avoidance_r</th>\n",
              "      <th>...</th>\n",
              "      <th>sex</th>\n",
              "      <th>Site</th>\n",
              "      <th>smoke</th>\n",
              "      <th>socialdiversity</th>\n",
              "      <th>socialembedded</th>\n",
              "      <th>socTherm</th>\n",
              "      <th>soliTherm</th>\n",
              "      <th>stress</th>\n",
              "      <th>site_idx</th>\n",
              "      <th>obs_id</th>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>Site</th>\n",
              "      <th>obs_id</th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th rowspan=\"10\" valign=\"top\">Kassel</th>\n",
              "      <th>0</th>\n",
              "      <td>1955.0</td>\n",
              "      <td>3.500000</td>\n",
              "      <td>0.323789</td>\n",
              "      <td>0.0</td>\n",
              "      <td>3.222222</td>\n",
              "      <td>3.666667</td>\n",
              "      <td>89.0</td>\n",
              "      <td>36.300</td>\n",
              "      <td>3.777778</td>\n",
              "      <td>1.061474</td>\n",
              "      <td>...</td>\n",
              "      <td>1.0</td>\n",
              "      <td>Kassel</td>\n",
              "      <td>2.0</td>\n",
              "      <td>6</td>\n",
              "      <td>3</td>\n",
              "      <td>3.0</td>\n",
              "      <td>3.125</td>\n",
              "      <td>30</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>1959.0</td>\n",
              "      <td>1.500000</td>\n",
              "      <td>-1.380456</td>\n",
              "      <td>0.0</td>\n",
              "      <td>4.444444</td>\n",
              "      <td>2.111111</td>\n",
              "      <td>NaN</td>\n",
              "      <td>36.600</td>\n",
              "      <td>2.722222</td>\n",
              "      <td>-0.124879</td>\n",
              "      <td>...</td>\n",
              "      <td>2.0</td>\n",
              "      <td>Kassel</td>\n",
              "      <td>2.0</td>\n",
              "      <td>7</td>\n",
              "      <td>4</td>\n",
              "      <td>2.4</td>\n",
              "      <td>3.125</td>\n",
              "      <td>30</td>\n",
              "      <td>0</td>\n",
              "      <td>1</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>1965.0</td>\n",
              "      <td>1.333333</td>\n",
              "      <td>-1.522476</td>\n",
              "      <td>0.0</td>\n",
              "      <td>3.666667</td>\n",
              "      <td>2.000000</td>\n",
              "      <td>NaN</td>\n",
              "      <td>35.450</td>\n",
              "      <td>2.611111</td>\n",
              "      <td>-0.249759</td>\n",
              "      <td>...</td>\n",
              "      <td>1.0</td>\n",
              "      <td>Kassel</td>\n",
              "      <td>2.0</td>\n",
              "      <td>5</td>\n",
              "      <td>2</td>\n",
              "      <td>3.8</td>\n",
              "      <td>2.375</td>\n",
              "      <td>31</td>\n",
              "      <td>0</td>\n",
              "      <td>2</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>1966.0</td>\n",
              "      <td>3.222222</td>\n",
              "      <td>0.087088</td>\n",
              "      <td>0.0</td>\n",
              "      <td>4.333333</td>\n",
              "      <td>1.000000</td>\n",
              "      <td>NaN</td>\n",
              "      <td>36.630</td>\n",
              "      <td>3.166667</td>\n",
              "      <td>0.374638</td>\n",
              "      <td>...</td>\n",
              "      <td>2.0</td>\n",
              "      <td>Kassel</td>\n",
              "      <td>2.0</td>\n",
              "      <td>9</td>\n",
              "      <td>2</td>\n",
              "      <td>3.6</td>\n",
              "      <td>2.875</td>\n",
              "      <td>47</td>\n",
              "      <td>0</td>\n",
              "      <td>3</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>1969.0</td>\n",
              "      <td>3.444444</td>\n",
              "      <td>0.276448</td>\n",
              "      <td>0.0</td>\n",
              "      <td>3.555556</td>\n",
              "      <td>1.444444</td>\n",
              "      <td>69.0</td>\n",
              "      <td>36.465</td>\n",
              "      <td>3.222222</td>\n",
              "      <td>0.437077</td>\n",
              "      <td>...</td>\n",
              "      <td>2.0</td>\n",
              "      <td>Kassel</td>\n",
              "      <td>2.0</td>\n",
              "      <td>11</td>\n",
              "      <td>3</td>\n",
              "      <td>2.2</td>\n",
              "      <td>2.750</td>\n",
              "      <td>50</td>\n",
              "      <td>0</td>\n",
              "      <td>4</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>5</th>\n",
              "      <td>1976.0</td>\n",
              "      <td>2.444444</td>\n",
              "      <td>-0.575674</td>\n",
              "      <td>0.0</td>\n",
              "      <td>4.555556</td>\n",
              "      <td>1.444444</td>\n",
              "      <td>73.0</td>\n",
              "      <td>36.500</td>\n",
              "      <td>2.222222</td>\n",
              "      <td>-0.686836</td>\n",
              "      <td>...</td>\n",
              "      <td>2.0</td>\n",
              "      <td>Kassel</td>\n",
              "      <td>2.0</td>\n",
              "      <td>12</td>\n",
              "      <td>3</td>\n",
              "      <td>3.2</td>\n",
              "      <td>2.625</td>\n",
              "      <td>42</td>\n",
              "      <td>0</td>\n",
              "      <td>5</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>6</th>\n",
              "      <td>1979.0</td>\n",
              "      <td>3.500000</td>\n",
              "      <td>0.323789</td>\n",
              "      <td>0.0</td>\n",
              "      <td>4.111111</td>\n",
              "      <td>2.777778</td>\n",
              "      <td>68.0</td>\n",
              "      <td>35.300</td>\n",
              "      <td>3.555556</td>\n",
              "      <td>0.811715</td>\n",
              "      <td>...</td>\n",
              "      <td>1.0</td>\n",
              "      <td>Kassel</td>\n",
              "      <td>2.0</td>\n",
              "      <td>7</td>\n",
              "      <td>2</td>\n",
              "      <td>2.8</td>\n",
              "      <td>4.125</td>\n",
              "      <td>37</td>\n",
              "      <td>0</td>\n",
              "      <td>6</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>7</th>\n",
              "      <td>1979.0</td>\n",
              "      <td>1.444444</td>\n",
              "      <td>-1.427796</td>\n",
              "      <td>0.0</td>\n",
              "      <td>4.777778</td>\n",
              "      <td>3.000000</td>\n",
              "      <td>54.0</td>\n",
              "      <td>36.100</td>\n",
              "      <td>2.222222</td>\n",
              "      <td>-0.686836</td>\n",
              "      <td>...</td>\n",
              "      <td>2.0</td>\n",
              "      <td>Kassel</td>\n",
              "      <td>2.0</td>\n",
              "      <td>8</td>\n",
              "      <td>2</td>\n",
              "      <td>3.4</td>\n",
              "      <td>2.750</td>\n",
              "      <td>33</td>\n",
              "      <td>0</td>\n",
              "      <td>7</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>8</th>\n",
              "      <td>1980.0</td>\n",
              "      <td>2.888889</td>\n",
              "      <td>-0.196953</td>\n",
              "      <td>0.0</td>\n",
              "      <td>4.111111</td>\n",
              "      <td>2.333333</td>\n",
              "      <td>NaN</td>\n",
              "      <td>36.400</td>\n",
              "      <td>3.388889</td>\n",
              "      <td>0.624396</td>\n",
              "      <td>...</td>\n",
              "      <td>1.0</td>\n",
              "      <td>Kassel</td>\n",
              "      <td>1.0</td>\n",
              "      <td>9</td>\n",
              "      <td>4</td>\n",
              "      <td>2.8</td>\n",
              "      <td>3.750</td>\n",
              "      <td>42</td>\n",
              "      <td>0</td>\n",
              "      <td>8</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>9</th>\n",
              "      <td>1981.0</td>\n",
              "      <td>4.555556</td>\n",
              "      <td>1.223251</td>\n",
              "      <td>0.0</td>\n",
              "      <td>4.777778</td>\n",
              "      <td>3.333333</td>\n",
              "      <td>61.0</td>\n",
              "      <td>36.750</td>\n",
              "      <td>4.000000</td>\n",
              "      <td>1.311232</td>\n",
              "      <td>...</td>\n",
              "      <td>2.0</td>\n",
              "      <td>Kassel</td>\n",
              "      <td>2.0</td>\n",
              "      <td>9</td>\n",
              "      <td>3</td>\n",
              "      <td>2.6</td>\n",
              "      <td>3.875</td>\n",
              "      <td>40</td>\n",
              "      <td>0</td>\n",
              "      <td>9</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "<p>10 rows × 39 columns</p>\n",
              "</div>"
            ],
            "text/plain": [
              "                  age   anxiety  anxiety_r  artgluctot  attachhome  \\\n",
              "Site   obs_id                                                        \n",
              "Kassel 0       1955.0  3.500000   0.323789         0.0    3.222222   \n",
              "       1       1959.0  1.500000  -1.380456         0.0    4.444444   \n",
              "       2       1965.0  1.333333  -1.522476         0.0    3.666667   \n",
              "       3       1966.0  3.222222   0.087088         0.0    4.333333   \n",
              "       4       1969.0  3.444444   0.276448         0.0    3.555556   \n",
              "       5       1976.0  2.444444  -0.575674         0.0    4.555556   \n",
              "       6       1979.0  3.500000   0.323789         0.0    4.111111   \n",
              "       7       1979.0  1.444444  -1.427796         0.0    4.777778   \n",
              "       8       1980.0  2.888889  -0.196953         0.0    4.111111   \n",
              "       9       1981.0  4.555556   1.223251         0.0    4.777778   \n",
              "\n",
              "               attachphone  AvgHumidity  avgtemp  avoidance  avoidance_r  ...  \\\n",
              "Site   obs_id                                                             ...   \n",
              "Kassel 0          3.666667         89.0   36.300   3.777778     1.061474  ...   \n",
              "       1          2.111111          NaN   36.600   2.722222    -0.124879  ...   \n",
              "       2          2.000000          NaN   35.450   2.611111    -0.249759  ...   \n",
              "       3          1.000000          NaN   36.630   3.166667     0.374638  ...   \n",
              "       4          1.444444         69.0   36.465   3.222222     0.437077  ...   \n",
              "       5          1.444444         73.0   36.500   2.222222    -0.686836  ...   \n",
              "       6          2.777778         68.0   35.300   3.555556     0.811715  ...   \n",
              "       7          3.000000         54.0   36.100   2.222222    -0.686836  ...   \n",
              "       8          2.333333          NaN   36.400   3.388889     0.624396  ...   \n",
              "       9          3.333333         61.0   36.750   4.000000     1.311232  ...   \n",
              "\n",
              "               sex    Site  smoke  socialdiversity  socialembedded  socTherm  \\\n",
              "Site   obs_id                                                                  \n",
              "Kassel 0       1.0  Kassel    2.0                6               3       3.0   \n",
              "       1       2.0  Kassel    2.0                7               4       2.4   \n",
              "       2       1.0  Kassel    2.0                5               2       3.8   \n",
              "       3       2.0  Kassel    2.0                9               2       3.6   \n",
              "       4       2.0  Kassel    2.0               11               3       2.2   \n",
              "       5       2.0  Kassel    2.0               12               3       3.2   \n",
              "       6       1.0  Kassel    2.0                7               2       2.8   \n",
              "       7       2.0  Kassel    2.0                8               2       3.4   \n",
              "       8       1.0  Kassel    1.0                9               4       2.8   \n",
              "       9       2.0  Kassel    2.0                9               3       2.6   \n",
              "\n",
              "               soliTherm  stress  site_idx  obs_id  \n",
              "Site   obs_id                                       \n",
              "Kassel 0           3.125      30         0       0  \n",
              "       1           3.125      30         0       1  \n",
              "       2           2.375      31         0       2  \n",
              "       3           2.875      47         0       3  \n",
              "       4           2.750      50         0       4  \n",
              "       5           2.625      42         0       5  \n",
              "       6           4.125      37         0       6  \n",
              "       7           2.750      33         0       7  \n",
              "       8           3.750      42         0       8  \n",
              "       9           3.875      40         0       9  \n",
              "\n",
              "[10 rows x 39 columns]"
            ]
          },
          "execution_count": 4,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# 选取5个站点\n",
        "first5_site = ['Southampton','METU','Kassel','Tsinghua','Oslo']\n",
        "df_first5 = df_raw.query(\"Site in @first5_site\")\n",
        "\n",
        "#为站点生成索引，为被试生成索引\n",
        "df_first5[\"site_idx\"] = pd.factorize(df_first5.Site)[0]\n",
        "df_first5[\"obs_id\"] = range(len(df_first5))\n",
        "\n",
        "#设置索引，方便之后调用数据\n",
        "df_first5.set_index(['Site','obs_id'],inplace=True,drop=False)\n",
        "df_first5.head(10)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "2D86124910AE4390B74939FFA27A829C",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 数据的层级结构  \n",
        "\n",
        "![Image Name](https://cdn.kesci.com/upload/s5eeyrh0s5.png?imageView2/0/w/960/h/960)  \n",
        "\n",
        "\n",
        "* $j$来表示站点，$j \\in \\{1,2, \\ldots, 5\\}$  \n",
        "* $i$来表示站点内部的每一个数据$i \\in \\{1,2,\\ldots,n_j\\}$  \n",
        "* 每一个被试的数据可以被表示为$Y_{ij}$，表示站点$j$内的第$i$个被试的自我控制分数观测值  \n",
        "\n",
        "$$  \n",
        "Y := \\left((Y_{11}, Y_{21}, \\ldots, Y_{n_1,1}), (Y_{12}, Y_{22}, \\ldots, Y_{n_2,2}), \\ldots, (Y_{1,5}, Y_{2,5}, \\ldots, Y_{n_{5},5})\\right)  .  \n",
        "$$  \n",
        "\n",
        "🤔可以想象，如果我们忽略了数据的分组结构，则可能违反模型使用的前提。  \n",
        "\n",
        "* 例如，由于反应来自同一个被试，数据点之间是相互关联的，在对被试的反应时使用回归模型时，需要考虑独立性假设是否被违反。  "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "34812936D2EF479692B2B3226ABBEE58",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "面对这样的层级数据，我们有三种建模思路(对应三种不同的假设)：  \n",
        "- 完全池化(complete pooling)  \n",
        "- 非池化(no pooling)  \n",
        "- 部分池化(patial pooling)  "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "ADDED191216C409FBD62217A46F0344A",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## 完全池化(Complete pooling)模型  \n",
        "\n",
        "* 在完全池化模型中，我们忽略个体的分组信息(Site)，认为这些个体直接来自一个更大的总体  \n",
        "\n",
        "![Image Name](https://cdn.kesci.com/upload/s5eflzcpoj.png?imageView2/0/w/960/h/960)  \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "F53C6E1ABDA24961A5A7546E41701125",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "我们可以使用正态模型来建立 complete pooled 模型：  \n",
        "\n",
        "$$  \n",
        "\\begin{split}  \n",
        "Y_{ij} | \\mu, \\sigma & \\sim N(\\mu, \\sigma^2) \\\\  \n",
        "\\mu    & \\sim N(0, 50^2) \\\\  \n",
        "\\sigma & \\sim \\text{Exp}(1) \\\\  \n",
        "\\end{split}  \n",
        "$$  \n",
        "\n",
        "**global parameters $\\theta$**  \n",
        "\n",
        "* $\\mu$为自我控制分数在**总体**中的均值(global mean)  \n",
        "\n",
        "* $\\sigma$为自我控制分数均值在**总体**中的标准差(global standard deviation)  \n",
        "\n",
        "![Image Name](https://cdn.kesci.com/upload/image/rmcx1f4dlx.png?imageView2/0/w/960/h/960)  \n",
        "\n",
        "> 下图给出了所有552个被试自我控制分数的分布情况，分数范围大概在20-60之间"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "collapsed": false,
        "id": "DF9F41D471214AF0BD1FF30468A88A6E",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/187F763E48114CAA83EEF80339913368/s5jqbxw9rg.png\">"
            ],
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "data_scontrol = df_first5[\"scontrol\"]\n",
        "# 绘制直方图\n",
        "ax = sns.histplot(data_scontrol, kde=True, bins=30)\n",
        "\n",
        "# 计算并画出均值线\n",
        "mean = np.mean(data_scontrol)\n",
        "ax.axvline(mean, color='red', linestyle='--')\n",
        "\n",
        "# 在直方图上添加散点以标记每个样本\n",
        "ax = sns.rugplot(data_scontrol, color=\"orange\", height=0.1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "2D0EB354AB914D24BC9225509700D467",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 模型定义及MCMC采样  \n",
        "\n",
        "$$  \n",
        "\\begin{split}  \n",
        "Y_{ij} | \\mu, \\sigma & \\sim N(\\mu, \\sigma^2) \\\\  \n",
        "\\mu    & \\sim N(0, 50^2) \\\\  \n",
        "\\sigma & \\sim \\text{Exp}(1) \\\\  \n",
        "\\end{split}  \n",
        "$$  \n",
        "\n",
        "根据公式使用 pymc 定义模型："
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "collapsed": false,
        "id": "21B29FC65FB940A689606E0B5995CF7A",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Auto-assigning NUTS sampler...\n",
            "Initializing NUTS using jitter+adapt_diag...\n",
            "Multiprocess sampling (4 chains in 4 jobs)\n",
            "NUTS: [mu, sigma]\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 2 seconds.\n"
          ]
        }
      ],
      "source": [
        "with pm.Model() as complete_pooled_model:\n",
        "\n",
        "    #定义beta_0\n",
        "    mu = pm.Normal(\"mu\", mu=0, sigma=50)  \n",
        "    #定义sigma                  \n",
        "    sigma = pm.Exponential(\"sigma\", 1)           \n",
        "    #定义似然：预测值y符合N(mu, sigma)分布；传入实际数据y 自我控制水平 df_first5.scontrol\n",
        "    likelihood = pm.Normal(\"y_est\", mu=mu, sigma=sigma, observed=df_first5.scontrol)   \n",
        "    # 进行采样，默认为 chains=4, samples=1000,burn=1000\n",
        "    complete_trace = pm.sample(random_seed=84735)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "collapsed": false,
        "id": "A1C27E7C5E404C748EDE670CB37CAC6F",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/46D2FE630C1C4D3F874DC434ECE19D80/s5jqcm660v.svg\">"
            ],
            "text/plain": [
              "<graphviz.graphs.Digraph at 0x7f985dab9370>"
            ]
          },
          "execution_count": 7,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "pm.model_to_graphviz(complete_pooled_model)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "38783D73191B4BC98086F8328104EED4",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 查看后验参数估计"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "collapsed": false,
        "id": "9D2AB6F859F94FCB85F293ABAA896709",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[<Axes: title={'center': 'mu'}>, <Axes: title={'center': 'mu'}>],\n",
              "       [<Axes: title={'center': 'sigma'}>,\n",
              "        <Axes: title={'center': 'sigma'}>]], dtype=object)"
            ]
          },
          "execution_count": 8,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/A584AA2375C142279A67901611BDD4AE/s5jqcp46un.png\">"
            ],
            "text/plain": [
              "<Figure size 1500x600 with 4 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "az.plot_trace(complete_trace,\n",
        "              compact=False,\n",
        "              figsize=(15,6))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "collapsed": false,
        "id": "E5E46A8728974EACB7016F5F3E184991",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>mean</th>\n",
              "      <th>sd</th>\n",
              "      <th>hdi_3%</th>\n",
              "      <th>hdi_97%</th>\n",
              "      <th>mcse_mean</th>\n",
              "      <th>mcse_sd</th>\n",
              "      <th>ess_bulk</th>\n",
              "      <th>ess_tail</th>\n",
              "      <th>r_hat</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>mu</th>\n",
              "      <td>40.448</td>\n",
              "      <td>0.324</td>\n",
              "      <td>39.864</td>\n",
              "      <td>41.070</td>\n",
              "      <td>0.005</td>\n",
              "      <td>0.004</td>\n",
              "      <td>4071.0</td>\n",
              "      <td>2795.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>sigma</th>\n",
              "      <td>7.513</td>\n",
              "      <td>0.225</td>\n",
              "      <td>7.116</td>\n",
              "      <td>7.941</td>\n",
              "      <td>0.003</td>\n",
              "      <td>0.002</td>\n",
              "      <td>4163.0</td>\n",
              "      <td>2928.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "         mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  \\\n",
              "mu     40.448  0.324  39.864   41.070      0.005    0.004    4071.0    2795.0   \n",
              "sigma   7.513  0.225   7.116    7.941      0.003    0.002    4163.0    2928.0   \n",
              "\n",
              "       r_hat  \n",
              "mu       1.0  \n",
              "sigma    1.0  "
            ]
          },
          "execution_count": 9,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "az.summary(complete_trace)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "66AB3DC0FC2448BE995997CDAD4AF0BB",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 后验预测分布  \n",
        "\n",
        "* 我们可以画出所有预测值的95%后验预测可信区间  \n",
        "\n",
        "* 同时可以观察真实值落在可信区间的情况  \n",
        "\n",
        "> (在lec11中，我们介绍过可信区间的绘制，主要借助`az.summary`生成后验预测总结的表格，并判断真实值是否落在可信区间内)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "collapsed": false,
        "id": "AFDC79A159F8441791021984075BAB46",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Sampling: [y_est]\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# 进行后验预测\n",
        "complete_ppc = pm.sample_posterior_predictive(complete_trace,\n",
        "                                              model=complete_pooled_model)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "collapsed": false,
        "id": "C5E8386680A64F2EB28C4D9B3CA0050E",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>mean</th>\n",
              "      <th>sd</th>\n",
              "      <th>hdi_2.5%</th>\n",
              "      <th>hdi_97.5%</th>\n",
              "      <th>mcse_mean</th>\n",
              "      <th>mcse_sd</th>\n",
              "      <th>ess_bulk</th>\n",
              "      <th>ess_tail</th>\n",
              "      <th>r_hat</th>\n",
              "      <th>obs_id</th>\n",
              "      <th>y</th>\n",
              "      <th>site</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>y_est[0]</th>\n",
              "      <td>40.571</td>\n",
              "      <td>7.511</td>\n",
              "      <td>25.067</td>\n",
              "      <td>54.316</td>\n",
              "      <td>0.118</td>\n",
              "      <td>0.084</td>\n",
              "      <td>4023.0</td>\n",
              "      <td>3757.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>0</td>\n",
              "      <td>47</td>\n",
              "      <td>Kassel</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>y_est[1]</th>\n",
              "      <td>40.580</td>\n",
              "      <td>7.531</td>\n",
              "      <td>24.656</td>\n",
              "      <td>53.966</td>\n",
              "      <td>0.117</td>\n",
              "      <td>0.083</td>\n",
              "      <td>4118.0</td>\n",
              "      <td>3892.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>1</td>\n",
              "      <td>44</td>\n",
              "      <td>Kassel</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>y_est[2]</th>\n",
              "      <td>40.495</td>\n",
              "      <td>7.493</td>\n",
              "      <td>26.089</td>\n",
              "      <td>55.060</td>\n",
              "      <td>0.119</td>\n",
              "      <td>0.084</td>\n",
              "      <td>3975.0</td>\n",
              "      <td>3572.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>2</td>\n",
              "      <td>47</td>\n",
              "      <td>Kassel</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>y_est[3]</th>\n",
              "      <td>40.373</td>\n",
              "      <td>7.552</td>\n",
              "      <td>26.083</td>\n",
              "      <td>56.026</td>\n",
              "      <td>0.120</td>\n",
              "      <td>0.085</td>\n",
              "      <td>3997.0</td>\n",
              "      <td>3666.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>3</td>\n",
              "      <td>37</td>\n",
              "      <td>Kassel</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>y_est[4]</th>\n",
              "      <td>40.492</td>\n",
              "      <td>7.466</td>\n",
              "      <td>26.214</td>\n",
              "      <td>55.344</td>\n",
              "      <td>0.118</td>\n",
              "      <td>0.083</td>\n",
              "      <td>4020.0</td>\n",
              "      <td>3813.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>4</td>\n",
              "      <td>33</td>\n",
              "      <td>Kassel</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>...</th>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>y_est[547]</th>\n",
              "      <td>40.641</td>\n",
              "      <td>7.529</td>\n",
              "      <td>25.456</td>\n",
              "      <td>54.511</td>\n",
              "      <td>0.122</td>\n",
              "      <td>0.086</td>\n",
              "      <td>3835.0</td>\n",
              "      <td>3736.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>547</td>\n",
              "      <td>52</td>\n",
              "      <td>Tsinghua</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>y_est[548]</th>\n",
              "      <td>40.428</td>\n",
              "      <td>7.490</td>\n",
              "      <td>26.105</td>\n",
              "      <td>54.716</td>\n",
              "      <td>0.119</td>\n",
              "      <td>0.084</td>\n",
              "      <td>3946.0</td>\n",
              "      <td>3793.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>548</td>\n",
              "      <td>33</td>\n",
              "      <td>Tsinghua</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>y_est[549]</th>\n",
              "      <td>40.425</td>\n",
              "      <td>7.497</td>\n",
              "      <td>25.979</td>\n",
              "      <td>55.412</td>\n",
              "      <td>0.118</td>\n",
              "      <td>0.083</td>\n",
              "      <td>4035.0</td>\n",
              "      <td>3808.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>549</td>\n",
              "      <td>35</td>\n",
              "      <td>Tsinghua</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>y_est[550]</th>\n",
              "      <td>40.872</td>\n",
              "      <td>7.521</td>\n",
              "      <td>26.337</td>\n",
              "      <td>55.374</td>\n",
              "      <td>0.116</td>\n",
              "      <td>0.083</td>\n",
              "      <td>4173.0</td>\n",
              "      <td>4015.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>550</td>\n",
              "      <td>38</td>\n",
              "      <td>Tsinghua</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>y_est[551]</th>\n",
              "      <td>40.857</td>\n",
              "      <td>7.586</td>\n",
              "      <td>25.969</td>\n",
              "      <td>55.533</td>\n",
              "      <td>0.119</td>\n",
              "      <td>0.084</td>\n",
              "      <td>4046.0</td>\n",
              "      <td>3960.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>551</td>\n",
              "      <td>39</td>\n",
              "      <td>Tsinghua</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "<p>552 rows × 12 columns</p>\n",
              "</div>"
            ],
            "text/plain": [
              "              mean     sd  hdi_2.5%  hdi_97.5%  mcse_mean  mcse_sd  ess_bulk  \\\n",
              "y_est[0]    40.571  7.511    25.067     54.316      0.118    0.084    4023.0   \n",
              "y_est[1]    40.580  7.531    24.656     53.966      0.117    0.083    4118.0   \n",
              "y_est[2]    40.495  7.493    26.089     55.060      0.119    0.084    3975.0   \n",
              "y_est[3]    40.373  7.552    26.083     56.026      0.120    0.085    3997.0   \n",
              "y_est[4]    40.492  7.466    26.214     55.344      0.118    0.083    4020.0   \n",
              "...            ...    ...       ...        ...        ...      ...       ...   \n",
              "y_est[547]  40.641  7.529    25.456     54.511      0.122    0.086    3835.0   \n",
              "y_est[548]  40.428  7.490    26.105     54.716      0.119    0.084    3946.0   \n",
              "y_est[549]  40.425  7.497    25.979     55.412      0.118    0.083    4035.0   \n",
              "y_est[550]  40.872  7.521    26.337     55.374      0.116    0.083    4173.0   \n",
              "y_est[551]  40.857  7.586    25.969     55.533      0.119    0.084    4046.0   \n",
              "\n",
              "            ess_tail  r_hat  obs_id   y      site  \n",
              "y_est[0]      3757.0    1.0       0  47    Kassel  \n",
              "y_est[1]      3892.0    1.0       1  44    Kassel  \n",
              "y_est[2]      3572.0    1.0       2  47    Kassel  \n",
              "y_est[3]      3666.0    1.0       3  37    Kassel  \n",
              "y_est[4]      3813.0    1.0       4  33    Kassel  \n",
              "...              ...    ...     ...  ..       ...  \n",
              "y_est[547]    3736.0    1.0     547  52  Tsinghua  \n",
              "y_est[548]    3793.0    1.0     548  33  Tsinghua  \n",
              "y_est[549]    3808.0    1.0     549  35  Tsinghua  \n",
              "y_est[550]    4015.0    1.0     550  38  Tsinghua  \n",
              "y_est[551]    3960.0    1.0     551  39  Tsinghua  \n",
              "\n",
              "[552 rows x 12 columns]"
            ]
          },
          "execution_count": 11,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# 定义函数，计算 95%hdi\n",
        "def ppc_sum(ppc, data):\n",
        "    \n",
        "    hdi_sum = az.summary(ppc, hdi_prob=0.95)\n",
        "    hdi_sum[\"obs_id\"] = data[\"obs_id\"].values\n",
        "    hdi_sum[\"y\"] = data[\"scontrol\"].values\n",
        "    hdi_sum[\"site\"] = data[\"Site\"].values\n",
        "\n",
        "    return hdi_sum\n",
        "\n",
        "# 计算后验预测的 95%hdi\n",
        "complete_hdi_sum = ppc_sum(ppc = complete_ppc, data=df_first5)\n",
        "complete_hdi_sum"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {
        "collapsed": false,
        "id": "818F9910FA2A41D6A6FA3342FF031A26",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 定义函数绘制超出 95%hdi 的点\n",
        "from matplotlib.lines import Line2D\n",
        "\n",
        "def ppc_plot(hdi_sum):\n",
        "    fig, ax =  plt.subplots(figsize=(15,6))\n",
        "\n",
        "    #生成颜色条件，根据站点生成不同的颜色（可信区间）\n",
        "    unique_sites = hdi_sum[\"site\"].unique()\n",
        "    conditions=[]\n",
        "    colors=[]\n",
        "    for i, site in enumerate(unique_sites):\n",
        "        condition = hdi_sum[\"site\"] == site\n",
        "        conditions.append(condition)\n",
        "        color = f\"C{i}\"\n",
        "        colors.append(color)\n",
        "        \n",
        "    hdi_colors = np.select(conditions,colors)\n",
        "    #绘制94%的可信区间\n",
        "    HDI = ax.vlines(hdi_sum[\"obs_id\"], \n",
        "            hdi_sum[\"hdi_2.5%\"], hdi_sum[\"hdi_97.5%\"], \n",
        "            color=hdi_colors, \n",
        "            alpha=0.5,\n",
        "            label=\"94% HDI\")\n",
        "    #绘制后验预测均值\n",
        "    pos_mean = ax.scatter(hdi_sum[\"obs_id\"], hdi_sum[\"mean\"],\n",
        "            marker=\"_\",\n",
        "            c = 'black',\n",
        "            alpha=0.2,\n",
        "            zorder = 2,\n",
        "            label=\"Posterior mean\")\n",
        "    #根据是否落在可信区间内选择不同的颜色\n",
        "    colors = np.where((hdi_sum[\"y\"] >= hdi_sum[\"hdi_2.5%\"]) & (hdi_sum[\"y\"] <= hdi_sum[\"hdi_97.5%\"]), \n",
        "                    '#2F5597', '#C00000')\n",
        "    #绘制真实值\n",
        "    ax.scatter(hdi_sum[\"obs_id\"], hdi_sum[\"y\"],\n",
        "            c = colors,\n",
        "            alpha=0.7,\n",
        "            zorder = 2)\n",
        "    # 设置图例的颜色、形状、名称\n",
        "    legend_color = ['#2F5597', '#C00000']\n",
        "    handles = [plt.Line2D([0], [0], \n",
        "                        marker='o', \n",
        "                        color='w', \n",
        "                        markerfacecolor=color, markersize=10) for color in legend_color]\n",
        "    handles += [HDI]\n",
        "    handles += [pos_mean]\n",
        "    labels = ['Within HDI', 'Outside HDI','94%HDI','Posterior mean']\n",
        "\n",
        "    plt.legend(handles=handles, \n",
        "               labels=labels,\n",
        "               loc='upper right',\n",
        "               bbox_to_anchor=(1.08, 1))\n",
        "    # 设置x轴的刻度，根据每个类别的数量确定刻度位置\n",
        "    count_per_site = hdi_sum.groupby(\"site\").size().values\n",
        "    cumulative_count = count_per_site.cumsum()\n",
        "    xtick = cumulative_count - count_per_site / 2\n",
        "    plt.xticks(xtick, hdi_sum[\"site\"].unique())\n",
        "\n",
        "    sns.despine()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "CAE781578DEE4FAF98E5BA56957A2C0C",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "由于我们省略了分组信息，假测所有的观测值都来自同一正态分布，因此所有观测值的后验预测均值都是相似的  \n",
        "* 纵坐标为自我控制分数的范围  \n",
        "* 横坐标代表每个观测值的排序  \n",
        "* 横线部分表示每个观测值的后验预测均值"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "collapsed": false,
        "id": "402C7CE11592459DA7870D618B09D025",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/5F13D450B7E243FF92A4896607DBC279/s5jqctfznk.png\">"
            ],
            "text/plain": [
              "<Figure size 1500x600 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "ppc_plot(hdi_sum=complete_hdi_sum)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "A2B484849A70496584E460A300D4BFD6",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## No pooling 非池化模型  \n",
        "\n",
        "在complete pooled模型中，我们忽略了数据来自不同的站点这一事实  \n",
        "\n",
        "现在我们考虑另外一种情况，我们**假设五个站点分别来自不同的分布，对五个站点进行不同的分析**  \n",
        "* 注意：我们假定，各站点数据之间完全没有关联，不同站点之间彼此独立。  \n",
        "* 从统计上讲，假定各站点之间的参数(例如均值$\\mu$)没有关系，或者说是完全异质。  "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {
        "collapsed": false,
        "id": "79D5DC5E9C7B4B39ACC1A57B02EC779F",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/6B7FC17256B94A048A5D66DC12E3D149/s5jqctq6b1.png\">"
            ],
            "text/plain": [
              "<Figure size 1800x700 with 2 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "fig, ax = plt.subplots(1,2,figsize=(18,7))\n",
        "sns.boxplot(data=df_first5,\n",
        "            x=\"Site\",\n",
        "            y=\"scontrol\",\n",
        "            ax=ax[0])\n",
        "\n",
        "sns.kdeplot(data=df_first5,\n",
        "            x=\"scontrol\",\n",
        "            hue=\"Site\",\n",
        "             ax=ax[1])\n",
        "sns.despine()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "8DBDAB083D8F484F9E815600F3F8B339",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### Group-specific parameters  \n",
        "\n",
        "* 在完全池化模型中，我们使用了正态分布的参数来自总体层面；  \n",
        "\n",
        "* 在非池化模型中，我们认为正态分布的参数在组与组之间是不同的(group-specific) ---- 可以认为么个站点的数据对应一个独立的正态分布模型  \n",
        "\n",
        "![Image Name](https://cdn.kesci.com/upload/s5eh3sji4k.png?imageView2/0/w/960/h/960)  \n",
        "\n",
        "1. 使用 $\\mu_j$ 来表示每个站点的自我控制分数均值  \n",
        "\n",
        "    * 不同站点的$\\mu_j$不同  \n",
        "\n",
        "    * 同一站点内的个体服从以 $\\mu_j$ 为均值的正态分布模型  \n",
        "\n",
        "3. 同样，使用 $\\sigma_j$ 来表示每个站点内部自我控制分数的变异性  \n",
        "\n",
        "    * 和 $\\mu_j$ 类似，不同站点的 $\\sigma_j$ 不同  \n",
        "    * 同一站点内的个体服从以 $\\sigma_j$ 为标准差的正态分布模型  \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "A840A0795F0249FC9622D912BAB22706",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "则对于站点 $j$ 内的个体来说，自我控制分数满足：  \n",
        "$$  \n",
        "Y_{ij} | \\mu_j, \\sigma \\sim N(\\mu_j, \\sigma_j^2) \\\\  \n",
        "\n",
        "\\mu_j  \\sim N(0, 50^2) \\\\  \n",
        "\n",
        "\\sigma_j \\sim \\text{Exp}(1) \\\\  \n",
        "$$  \n",
        "\n",
        "![Image Name](https://cdn.kesci.com/upload/s5eh3sji4k.png?imageView2/0/w/960/h/960)  "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "EF1026040150434DB4F0F20DA33F59FF",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 模型定义及MCMC采样  \n",
        "\n",
        "根据公式定义模型：  \n",
        "\n",
        "$$  \n",
        "Y_{ij} | \\mu_j, \\sigma \\sim N(\\mu_j, \\sigma_j^2) \\\\  \n",
        "\n",
        "\\mu_j  \\sim N(0, 50^2) \\\\  \n",
        "\n",
        "\\sigma_j \\sim \\text{Exp}(1) \\\\  \n",
        "$$  \n",
        "\n",
        "- 考虑到数据有5个站点，即 j = 1,2,3,4,5。因此，$\\mu_j$ 和 $\\sigma_j$ 也有5个值。  \n",
        "- 在pymc中，我们可以通过定义坐标 coords 来实现 `pm.Normal(..., dims=\"site\")  \n",
        "- 此外，每个个体的数据来自于某个站点，因此可以通过 `pm.MutableData(\"site\", df_first5.site_idx, dims=\"obs_id\")` 来定义个体数据 `obs_id` 和站点 `site` 之间的映射。"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {
        "collapsed": false,
        "id": "9D458CE83BA643338D5C1AB58C5FCC7D",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Auto-assigning NUTS sampler...\n",
            "Initializing NUTS using jitter+adapt_diag...\n",
            "Multiprocess sampling (4 chains in 4 jobs)\n",
            "NUTS: [mu, sigma]\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 5 seconds.\n"
          ]
        }
      ],
      "source": [
        "coords = {\"site\": df_first5[\"Site\"].unique(),\n",
        "          \"obs_id\": df_first5.obs_id}\n",
        "\n",
        "with pm.Model(coords=coords) as no_pooled_model:\n",
        "\n",
        "    #定义mu，指定dims=\"site\"，生成不同的mu \n",
        "    mu = pm.Normal(\"mu\", mu=0, sigma=50, dims=\"site\")                  \n",
        "    #定义sigma，指定dims=\"site\"，生成不同的sigma\n",
        "    sigma = pm.Exponential(\"sigma\", 2, dims=\"site\")            \n",
        "    #获得观测值对应的站点映射\n",
        "    site = pm.MutableData(\"site\", df_first5.site_idx, dims=\"obs_id\") \n",
        "    # 定义 likelihood\n",
        "    likelihood = pm.Normal(\"y_est\", mu=mu[site], sigma=sigma[site], observed=df_first5.scontrol, dims=\"obs_id\")\n",
        "\n",
        "    no_trace = pm.sample(random_seed=84735)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "metadata": {
        "collapsed": false,
        "id": "3FCF3FC4C83D491796EA2A41C747F6F3",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/FD9B3266776B4B81BEE741503EF1229F/s5jqdfkqhz.svg\">"
            ],
            "text/plain": [
              "<graphviz.graphs.Digraph at 0x7f9859828ee0>"
            ]
          },
          "execution_count": 16,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "pm.model_to_graphviz(no_pooled_model)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "082EE78AF8AE458BBE4159449D3E6C53",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 查看后验参数估计  \n",
        "\n",
        "- 可以发现，对于每个站点，均有不同的参数值 (包括 $\\mu$ 和 $\\sigma$)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "metadata": {
        "collapsed": false,
        "id": "8FD3D43BAE634F849FB742D3DDEA7F5F",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/F4920CF066874D2C9CA97BE5CC410316/s5jqdsboct.png\">"
            ],
            "text/plain": [
              "<Figure size 2000x5000 with 20 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "ax = az.plot_trace(\n",
        "    no_trace,\n",
        "    compact=False,\n",
        "    figsize=(20,50))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "metadata": {
        "collapsed": false,
        "id": "066A298603DA49DCBC5B6D2015638AF7",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>mean</th>\n",
              "      <th>sd</th>\n",
              "      <th>hdi_3%</th>\n",
              "      <th>hdi_97%</th>\n",
              "      <th>mcse_mean</th>\n",
              "      <th>mcse_sd</th>\n",
              "      <th>ess_bulk</th>\n",
              "      <th>ess_tail</th>\n",
              "      <th>r_hat</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>mu[Kassel]</th>\n",
              "      <td>41.342</td>\n",
              "      <td>0.732</td>\n",
              "      <td>39.927</td>\n",
              "      <td>42.662</td>\n",
              "      <td>0.011</td>\n",
              "      <td>0.008</td>\n",
              "      <td>4710.0</td>\n",
              "      <td>2942.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>mu[METU]</th>\n",
              "      <td>39.854</td>\n",
              "      <td>0.602</td>\n",
              "      <td>38.768</td>\n",
              "      <td>41.016</td>\n",
              "      <td>0.009</td>\n",
              "      <td>0.006</td>\n",
              "      <td>4933.0</td>\n",
              "      <td>3081.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>mu[Oslo]</th>\n",
              "      <td>42.699</td>\n",
              "      <td>0.748</td>\n",
              "      <td>41.249</td>\n",
              "      <td>44.059</td>\n",
              "      <td>0.010</td>\n",
              "      <td>0.007</td>\n",
              "      <td>5493.0</td>\n",
              "      <td>3140.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>mu[Southampton]</th>\n",
              "      <td>38.592</td>\n",
              "      <td>1.900</td>\n",
              "      <td>34.871</td>\n",
              "      <td>42.057</td>\n",
              "      <td>0.027</td>\n",
              "      <td>0.019</td>\n",
              "      <td>4801.0</td>\n",
              "      <td>2991.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>mu[Tsinghua]</th>\n",
              "      <td>39.474</td>\n",
              "      <td>0.449</td>\n",
              "      <td>38.639</td>\n",
              "      <td>40.288</td>\n",
              "      <td>0.006</td>\n",
              "      <td>0.004</td>\n",
              "      <td>6037.0</td>\n",
              "      <td>3327.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>sigma[Kassel]</th>\n",
              "      <td>7.498</td>\n",
              "      <td>0.482</td>\n",
              "      <td>6.602</td>\n",
              "      <td>8.384</td>\n",
              "      <td>0.006</td>\n",
              "      <td>0.005</td>\n",
              "      <td>5740.0</td>\n",
              "      <td>3064.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>sigma[METU]</th>\n",
              "      <td>8.181</td>\n",
              "      <td>0.409</td>\n",
              "      <td>7.436</td>\n",
              "      <td>8.964</td>\n",
              "      <td>0.006</td>\n",
              "      <td>0.004</td>\n",
              "      <td>5387.0</td>\n",
              "      <td>3130.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>sigma[Oslo]</th>\n",
              "      <td>6.741</td>\n",
              "      <td>0.470</td>\n",
              "      <td>5.864</td>\n",
              "      <td>7.624</td>\n",
              "      <td>0.006</td>\n",
              "      <td>0.005</td>\n",
              "      <td>5644.0</td>\n",
              "      <td>3157.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>sigma[Southampton]</th>\n",
              "      <td>4.646</td>\n",
              "      <td>0.810</td>\n",
              "      <td>3.289</td>\n",
              "      <td>6.154</td>\n",
              "      <td>0.011</td>\n",
              "      <td>0.008</td>\n",
              "      <td>5578.0</td>\n",
              "      <td>2678.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>sigma[Tsinghua]</th>\n",
              "      <td>5.944</td>\n",
              "      <td>0.306</td>\n",
              "      <td>5.361</td>\n",
              "      <td>6.507</td>\n",
              "      <td>0.004</td>\n",
              "      <td>0.003</td>\n",
              "      <td>5130.0</td>\n",
              "      <td>2958.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "                      mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  \\\n",
              "mu[Kassel]          41.342  0.732  39.927   42.662      0.011    0.008   \n",
              "mu[METU]            39.854  0.602  38.768   41.016      0.009    0.006   \n",
              "mu[Oslo]            42.699  0.748  41.249   44.059      0.010    0.007   \n",
              "mu[Southampton]     38.592  1.900  34.871   42.057      0.027    0.019   \n",
              "mu[Tsinghua]        39.474  0.449  38.639   40.288      0.006    0.004   \n",
              "sigma[Kassel]        7.498  0.482   6.602    8.384      0.006    0.005   \n",
              "sigma[METU]          8.181  0.409   7.436    8.964      0.006    0.004   \n",
              "sigma[Oslo]          6.741  0.470   5.864    7.624      0.006    0.005   \n",
              "sigma[Southampton]   4.646  0.810   3.289    6.154      0.011    0.008   \n",
              "sigma[Tsinghua]      5.944  0.306   5.361    6.507      0.004    0.003   \n",
              "\n",
              "                    ess_bulk  ess_tail  r_hat  \n",
              "mu[Kassel]            4710.0    2942.0    1.0  \n",
              "mu[METU]              4933.0    3081.0    1.0  \n",
              "mu[Oslo]              5493.0    3140.0    1.0  \n",
              "mu[Southampton]       4801.0    2991.0    1.0  \n",
              "mu[Tsinghua]          6037.0    3327.0    1.0  \n",
              "sigma[Kassel]         5740.0    3064.0    1.0  \n",
              "sigma[METU]           5387.0    3130.0    1.0  \n",
              "sigma[Oslo]           5644.0    3157.0    1.0  \n",
              "sigma[Southampton]    5578.0    2678.0    1.0  \n",
              "sigma[Tsinghua]       5130.0    2958.0    1.0  "
            ]
          },
          "execution_count": 18,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "az.summary(no_trace)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "0E292DBB070946F097389839B40FC72C",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 后验预测分布"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "metadata": {
        "collapsed": false,
        "id": "3E4D80F435384EBE99451DB4C31BF78E",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Sampling: [y_est]\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "no_ppc = pm.sample_posterior_predictive(no_trace,\n",
        "                                        model=no_pooled_model)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "metadata": {
        "collapsed": false,
        "id": "5E7436DD611E458799F83482EAC35710",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "no_hdi_sum = ppc_sum(ppc = no_ppc,\n",
        "                data=df_first5)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "2AEE21BD9AF64F3C9F0BEBB26231AD21",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 非池化模型的缺点  \n",
        "\n",
        "可以看到在非池化模型中，每个组的均值与方差都是不同的，非池化模型充分考虑了每个组内部的情况，然而，这种模型的缺点可能包括：  \n",
        "\n",
        "1. 在小样本数据上，非池化模型存在过拟合的风险(如对于站点Southampton)  \n",
        "\n",
        "2. 非池化模型假设每个组都属于不同的分布，因此其得出的结果难以用来预测新组别的情况"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "metadata": {
        "collapsed": false,
        "id": "B15FEDCD58944348A807F1D7250FF8AA",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/BC293919183D45399FA71B86170B72AB/s5jqdtqx1t.png\">"
            ],
            "text/plain": [
              "<Figure size 1500x600 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "ppc_plot(hdi_sum=no_hdi_sum)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "98400C77FCEB4BC78E1F6C5AC17AD3DA",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## Partial pooling  \n",
        "\n",
        "- 完全池化模型仅考虑了个体差异，没有考虑站点之间的差异。  \n",
        "- 非池化模型考虑了站点的差异，但是容易受到站点内部数据的影响，忽视了站点间可能存在的关系。  \n",
        "\n",
        "部分池化 (Partial pooling)方法，是构建分层模型的关键。  \n",
        "- 它假设，不同站点 (layer2) 来自于一个关于站点的总体 (layer1)，站点形式的分布提供了对于组间变异(between variability)的解释。  \n",
        "- 此外，不同个体 (layer3) 又来自于不同站点 (layer2) ，每个站点内部个体形成的分布提供了对于组内变异(within variability)的解释。  \n",
        "\n",
        "![Image Name](https://cdn.kesci.com/upload/s5enwlkkz8.png?imageView2/0/w/960/h/960)  "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "E623501016214E12AE5F2BB2D836F519",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 对层级模型的定义  \n",
        "\n",
        "层级模型的数学形式：  \n",
        "\n",
        "$$  \n",
        "\\begin{array}{lrl}  \n",
        "\\text{Layer 1:} \\hspace{0.5in} Y_{ij} | \\mu_j, \\sigma_y \\sim N(\\mu_j, \\sigma_y^2) \\hspace{0.05in} \\text{根据站点均值生成个体数据}\\\\  \n",
        "\n",
        "\\text{Layer 2:} \\hspace{0.5in} \\mu_j | \\mu, \\sigma_\\mu \\stackrel{ind}{\\sim} N(\\mu, \\sigma_\\mu^2) \\hspace{0.13in} \\text{生成不同站点的均值}\\\\  \n",
        "\n",
        "\\text{Layer 3:}  \\hspace{0.5in}\\mu,\\sigma_y,\\sigma_\\mu \\hspace{0.85in}{超参数}\\\\  \n",
        "\\hspace{1in}\\mu \\sim N(40, 20^2) \\\\  \n",
        "\\hspace{1in}\\sigma_y \\sim \\text{Exp}(1) \\text{自我控制分数在组内的变异性} \\\\  \n",
        "\\hspace{1in} \\sigma_\\mu \\sim \\text{Exp}(1) \\text{均值在组间的变异性} \\\\  \n",
        "\\end{array}  \n",
        "$$  \n",
        "\n",
        "\n",
        "1. 在最上层：  \n",
        "    * 使用$\\mu$来表示总体的$Y$(global average)，即总体而言，自我控制分数的均值  \n",
        "    * 使用$\\sigma_\\mu$来表示组与组之间在$Y$均值上的变异性  \n",
        "    * 使用$\\sigma_y$来表示每个组内部$Y$的变异性(这里我们假设每个组内部的变异性相同)  \n",
        "2. 在第二层：  \n",
        "    * 使用$\\mu_j$来表示每个组$Y$的均值，而$\\mu_j$服从$N(\\mu, \\sigma_\\mu^2)$，这两个参数是从上一层得到的  \n",
        "3. 在第三层：  \n",
        "    * 使用$Y_{ij}$来表示组内每个个体的$Y$，而$Y_{ij}$服从$N(\\mu_j, \\sigma_y^2)$，这两个参数是从上一层得到的  \n",
        "\n",
        "![Image Name](https://cdn.kesci.com/upload/s5enwlkkz8.png?imageView2/0/w/500/h/500)  \n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "77958F236E614D24AD446D3DED9F23D1",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 另一种定义方式  \n",
        "\n",
        "* 我们也可以这样来理解层级模型：  \n",
        "\n",
        "    * $\\mu_j$表示：每个站点的自我控制均分  \n",
        "\n",
        "    * $\\sigma_\\mu$表示：$\\mu_j$偏离总体自我控制均分$\\mu$的程度  \n",
        "\n",
        "* 那么，$\\mu_j$与$\\mu$的关系式可以写为：  \n",
        "\n",
        "$$  \n",
        "\\mu_j = \\mu + b_{j}  \n",
        "$$  \n",
        "\n",
        "* 同时，$b_j$满足  \n",
        "$$  \n",
        "b_j \\sim N(0, \\sigma_\\mu^2)  \n",
        "$$  \n",
        "\n",
        "* 那么层级模型的Layer1 和 Layer2也可以写为：  \n",
        "\n",
        "$$  \n",
        "\\begin{split}  \n",
        "Y_{ij} | \\mu_j, \\sigma_y & \\sim N(\\mu_j, \\sigma_y^2) \\;\\; \\text{ with } \\;\\; \\mu_j = \\mu + b_{j}  \\\\  \n",
        "b_{j} | \\sigma_\\mu    & \\stackrel{ind}{\\sim} N(0, \\sigma_\\mu^2) \\\\  \n",
        "\\mu           & \\sim N(40, 20^2) \\\\  \n",
        "\\sigma_y      & \\sim \\text{Exp}(1) \\\\  \n",
        "\\sigma_\\mu    & \\sim \\text{Exp}(1) \\\\  \n",
        "\\end{split}  \n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "C30E8E4B6A2F4141A27E1CB401921A74",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 模型定义及MCMC采样  \n",
        "\n",
        "根据公式定义模型：  \n",
        "\n",
        "$$  \n",
        "\\begin{array}{lrl}  \n",
        "\\text{Layer 1:} \\hspace{0.5in} Y_{ij} | \\mu_j, \\sigma_y \\sim N(\\mu_j, \\sigma_y^2) \\hspace{0.05in} \\text{根据站点均值生成个体数据}\\\\  \n",
        "\n",
        "\\text{Layer 2:} \\hspace{0.5in} \\mu_j | \\mu, \\sigma_\\mu \\stackrel{ind}{\\sim} N(\\mu, \\sigma_\\mu^2) \\hspace{0.13in} \\text{生成不同站点的均值}\\\\  \n",
        "\n",
        "\\text{Layer 3:}  \\hspace{0.5in}\\mu,\\sigma_y,\\sigma_\\mu \\hspace{0.85in}{超参数}\\\\  \n",
        "\\hspace{1in}\\mu \\sim N(40, 20^2) \\\\  \n",
        "\\hspace{1in}\\sigma_y \\sim \\text{Exp}(1) \\text{自我控制分数在组内的变异性} \\\\  \n",
        "\\hspace{1in} \\sigma_\\mu \\sim \\text{Exp}(1) \\text{均值在组间的变异性} \\\\  \n",
        "\\end{array}  \n",
        "$$  \n",
        "\n",
        "- 相对于非池化模型。部分池化模型的关键在于定义，站点参数在总体上的变异，即 `mu = pm.Normal(\"mu\", mu=hyper_mu, sigma=var_mu, dims=\"site\")`  \n",
        "- 其次，是相似于非池化模型，个体数据在不同站点分布上的变异 `pm.Normal(\"y_est\", mu=mu[site], sigma=var_y, observed=df_first5.scontrol, dims=\"obs_id\")`  \n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "metadata": {
        "collapsed": false,
        "id": "15BAF93114DE4C8486A12187C5BB331B",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Auto-assigning NUTS sampler...\n",
            "Initializing NUTS using jitter+adapt_diag...\n",
            "Multiprocess sampling (4 chains in 4 jobs)\n",
            "NUTS: [within_variability, between_variability, hyper_mu, mu]\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 16 seconds.\n",
            "There were 414 divergences after tuning. Increase `target_accept` or reparameterize.\n"
          ]
        }
      ],
      "source": [
        "with pm.Model(coords=coords) as partial_pooled_model:\n",
        "    # Hyperpriors,定义全局参数\n",
        "    var_y = pm.Exponential(\"within_variability\", 1)\n",
        "    var_mu = pm.Exponential(\"between_variability\", 1)\n",
        "    hyper_mu = pm.Normal(\"hyper_mu\", mu=40, sigma=20)\n",
        "    # 定义站点参数\n",
        "    mu = pm.Normal(\"mu\", mu=hyper_mu, sigma=var_mu, dims=\"site\")\n",
        "    #获得观测值对应的站点映射\n",
        "    site = pm.MutableData(\"site\", df_first5.site_idx, dims=\"obs_id\")\n",
        "    # 定义 likelihood\n",
        "    likelihood = pm.Normal(\"y_est\", mu=mu[site], sigma=var_y, observed=df_first5.scontrol, dims=\"obs_id\")\n",
        "\n",
        "    partial_trace = pm.sample(draws=5000,                   # 使用mcmc方法进行采样，draws为采样次数\n",
        "                                tune=1000,                    # tune为调整采样策略的次数，可以决定这些结果是否要被保留\n",
        "                                chains=4,                     # 链数\n",
        "                                discard_tuned_samples= True,  # tune的结果将在采样结束后被丢弃\n",
        "                                random_seed=84735)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "metadata": {
        "collapsed": false,
        "id": "E56AC1FB0F564A0CB1953E5A345BE2FD",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/26C2FB8D34394BD394203357B4F8BE71/s5jqemcukh.svg\">"
            ],
            "text/plain": [
              "<graphviz.graphs.Digraph at 0x7f985df95e20>"
            ]
          },
          "execution_count": 23,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "pm.model_to_graphviz(partial_pooled_model)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {
        "collapsed": false,
        "id": "FA11BFB2B5454E8BA4440D88108941FC",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "with pm.Model(coords=coords) as partial_pooled_model_2:\n",
        "    # Hyperpriors\n",
        "    var_y = pm.Exponential(\"within_variability\", 1)\n",
        "    var_mu = pm.Exponential(\"between_variability\", 1)\n",
        "    hyper_mu = pm.Normal(\"hyper_mu\", mu=40, sigma=20)\n",
        "    mu_j = pm.Normal(\"mu_j\", mu=0, sigma=var_mu, dims=\"site\")\n",
        "    mu = pm.Deterministic(\"mu\", mu_j + hyper_mu)\n",
        "\n",
        "    site = pm.MutableData(\"site\", df_first5.site_idx, dims=\"obs_id\")\n",
        "\n",
        "    likelihood = pm.Normal(\"y_est\", mu=mu[site], sigma=var_y, observed=df_first5.scontrol, dims=\"obs_id\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "metadata": {
        "collapsed": false,
        "id": "02CB9483FC00430DBB6EA37C39258786",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/FA055610E44C43EB9EC08B284ED44DA0/s5jqempfpr.svg\">"
            ],
            "text/plain": [
              "<graphviz.graphs.Digraph at 0x7f985db686a0>"
            ]
          },
          "execution_count": 25,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "pm.model_to_graphviz(partial_pooled_model_2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "EC2E8A6611CA4B7B9047302103FA048B",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 查看后验参数估计  \n",
        "\n",
        "可以发现，  \n",
        "* 不同站点的均值 $\\mu[?]$ 接近站点的总体均值 $\\text{hyper}\\mu$  \n",
        "* 此外，组间变异和组内变异的结果也直观的显示了"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "metadata": {
        "collapsed": false,
        "id": "4CDBBACF0CF347CC9EB64EF06138A3F1",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>mean</th>\n",
              "      <th>sd</th>\n",
              "      <th>hdi_3%</th>\n",
              "      <th>hdi_97%</th>\n",
              "      <th>mcse_mean</th>\n",
              "      <th>mcse_sd</th>\n",
              "      <th>ess_bulk</th>\n",
              "      <th>ess_tail</th>\n",
              "      <th>r_hat</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>hyper_mu</th>\n",
              "      <td>40.629</td>\n",
              "      <td>0.729</td>\n",
              "      <td>39.337</td>\n",
              "      <td>42.100</td>\n",
              "      <td>0.007</td>\n",
              "      <td>0.005</td>\n",
              "      <td>10882.0</td>\n",
              "      <td>10954.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>mu[Kassel]</th>\n",
              "      <td>41.084</td>\n",
              "      <td>0.634</td>\n",
              "      <td>39.939</td>\n",
              "      <td>42.286</td>\n",
              "      <td>0.008</td>\n",
              "      <td>0.005</td>\n",
              "      <td>6816.0</td>\n",
              "      <td>8890.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>mu[METU]</th>\n",
              "      <td>40.055</td>\n",
              "      <td>0.513</td>\n",
              "      <td>39.070</td>\n",
              "      <td>41.009</td>\n",
              "      <td>0.005</td>\n",
              "      <td>0.004</td>\n",
              "      <td>9144.0</td>\n",
              "      <td>5299.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>mu[Oslo]</th>\n",
              "      <td>41.894</td>\n",
              "      <td>0.811</td>\n",
              "      <td>40.410</td>\n",
              "      <td>43.433</td>\n",
              "      <td>0.017</td>\n",
              "      <td>0.012</td>\n",
              "      <td>2231.0</td>\n",
              "      <td>742.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>mu[Southampton]</th>\n",
              "      <td>40.346</td>\n",
              "      <td>1.306</td>\n",
              "      <td>37.681</td>\n",
              "      <td>42.739</td>\n",
              "      <td>0.012</td>\n",
              "      <td>0.009</td>\n",
              "      <td>12110.0</td>\n",
              "      <td>10346.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>mu[Tsinghua]</th>\n",
              "      <td>39.773</td>\n",
              "      <td>0.544</td>\n",
              "      <td>38.722</td>\n",
              "      <td>40.757</td>\n",
              "      <td>0.005</td>\n",
              "      <td>0.004</td>\n",
              "      <td>10176.0</td>\n",
              "      <td>7878.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>within_variability</th>\n",
              "      <td>7.448</td>\n",
              "      <td>0.224</td>\n",
              "      <td>7.016</td>\n",
              "      <td>7.856</td>\n",
              "      <td>0.002</td>\n",
              "      <td>0.001</td>\n",
              "      <td>13778.0</td>\n",
              "      <td>11593.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>between_variability</th>\n",
              "      <td>1.211</td>\n",
              "      <td>0.635</td>\n",
              "      <td>0.237</td>\n",
              "      <td>2.344</td>\n",
              "      <td>0.012</td>\n",
              "      <td>0.008</td>\n",
              "      <td>1429.0</td>\n",
              "      <td>461.0</td>\n",
              "      <td>1.01</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "                       mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  \\\n",
              "hyper_mu             40.629  0.729  39.337   42.100      0.007    0.005   \n",
              "mu[Kassel]           41.084  0.634  39.939   42.286      0.008    0.005   \n",
              "mu[METU]             40.055  0.513  39.070   41.009      0.005    0.004   \n",
              "mu[Oslo]             41.894  0.811  40.410   43.433      0.017    0.012   \n",
              "mu[Southampton]      40.346  1.306  37.681   42.739      0.012    0.009   \n",
              "mu[Tsinghua]         39.773  0.544  38.722   40.757      0.005    0.004   \n",
              "within_variability    7.448  0.224   7.016    7.856      0.002    0.001   \n",
              "between_variability   1.211  0.635   0.237    2.344      0.012    0.008   \n",
              "\n",
              "                     ess_bulk  ess_tail  r_hat  \n",
              "hyper_mu              10882.0   10954.0   1.00  \n",
              "mu[Kassel]             6816.0    8890.0   1.00  \n",
              "mu[METU]               9144.0    5299.0   1.00  \n",
              "mu[Oslo]               2231.0     742.0   1.00  \n",
              "mu[Southampton]       12110.0   10346.0   1.00  \n",
              "mu[Tsinghua]          10176.0    7878.0   1.00  \n",
              "within_variability    13778.0   11593.0   1.00  \n",
              "between_variability    1429.0     461.0   1.01  "
            ]
          },
          "execution_count": 26,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "az.summary(partial_trace)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "metadata": {
        "collapsed": false,
        "id": "E61B431FBAA3468CAFF56647ADAD5FC7",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/C2E9CDD9357E4DE59B0061A11237AD36/s5jqexbo7a.png\">"
            ],
            "text/plain": [
              "<Figure size 2000x4000 with 16 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "with partial_pooled_model:\n",
        "    az.plot_trace(partial_trace,\n",
        "                  compact=False,\n",
        "                  figsize=(20,40))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "D02F629451D047C0952764CE846E4F0F",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## Within- vs between-group variability  \n",
        "\n",
        "* 在完全池化模型中，变异性的来源只有一种，即个体观测值 $Y$ 在总体中的变异  ($\\sigma$)  \n",
        "* 在非池化模型中，变异来源分散在不同站点中，即个体观测值在不同站点中的变异 ($\\sigma_j$)  \n",
        "\n",
        "在部分池化模型中，我们可以把个体观测值$Y$的变异来源分解成两个部分：  \n",
        "* 组内变异 $\\sigma_y^2$  \n",
        "* 组间变异 $\\sigma_{\\mu}^2$  \n",
        "\n",
        "$$  \n",
        "\\text{Var}(Y_{ij}) = \\sigma_y^2 + \\sigma_{\\mu}^2  \n",
        "$$  \n",
        "\n",
        "\n",
        "![Image Name](https://cdn.kesci.com/upload/s5eppwep1v.png?imageView2/0/w/600/h/600)  "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "518B47232AA7476EB871554F5B1A56F8",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "类似于单因素方差分析  \n",
        "\n",
        "$$  \n",
        "\\begin{split}  \n",
        "\\frac{\\sigma^2_y}{\\sigma^2_\\mu + \\sigma^2_y}  \n",
        "& = \\text{ $\\text{Var}(Y_{ij})$ 可以被组内方差解释的部分} \\\\  \n",
        "\\frac{\\sigma^2_\\mu}{\\sigma^2_\\mu + \\sigma^2_y}  \n",
        "& = \\text{$\\text{Var}(Y_{ij})$ 可以被组间方差解释的部分} \\\\  \n",
        "\\end{split}  \n",
        "$$  \n",
        "\n",
        "\n",
        "* 组间方差远大于组内方差($\\sigma_\\mu > \\sigma_y$)，那么组间方差可以解释大部分观测值的变异  \n",
        "* 组间方差远小于组内方差($\\sigma_\\mu < \\sigma_y$)，那么组与组之间的区别就不是很明显"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "7FC166111A104D4FA0C0CAFCE14A40F9",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "\n",
        "* 此外，若组间方差远大于组内方差，则可以说明组内的变异性很小，组内分数高度相关  \n",
        "\n",
        "$$  \n",
        "\\text{Cor}(Y_{ij}, Y_{kj}) = \\frac{\\sigma^2_\\mu}{\\sigma^2_\\mu + \\sigma^2_y}  \n",
        "$$  \n",
        "\n",
        "* 下图展示了三种组间方差与组内方差的分布情况(横轴为方差的大小，蓝色为组间方差，黑色为组内方差)  \n",
        "\n",
        "\n",
        "![Image Name](https://cdn.kesci.com/upload/s5k6q39qp6.png?imageView2/0/w/960/h/960)  \n",
        "\n",
        "\n",
        "- 组别越独特，σμ 相对越大，每**组内的相关性就越大**。  \n",
        "- 图 (a) 中，组间变异是组内变异的4倍，此时组内相关性为 0.8，接近于 1。  \n",
        "- 图 (c) 中，组间变异是组内变异的1/4倍，组内相关性为 0.2，接近 0。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "4C75DD46C668494E9D538CA2C994273F",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "我们可以使用`az.summary`来总结后验参数估计的情况，并计算组间组内变异和相关性"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 28,
      "metadata": {
        "collapsed": false,
        "id": "51563ECF911F4E13B08D10969ACA92A3",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "被组间方差所解释的部分： 0.025755900260321422\n",
            "被组内方差所解释的部分： 0.9742440997396786\n",
            "组内相关： 0.025755900260321422\n"
          ]
        }
      ],
      "source": [
        "# 提取组间和组内变异\n",
        "partial_para_sum = az.summary(partial_trace)\n",
        "between_sd = partial_para_sum.loc[\"between_variability\",\"mean\"]\n",
        "within_sd = partial_para_sum.loc[\"within_variability\",\"mean\"]\n",
        "# 计算变异占比\n",
        "var = between_sd**2 + within_sd**2\n",
        "print(\"被组间方差所解释的部分：\", between_sd**2/var)\n",
        "print(\"被组内方差所解释的部分：\", within_sd**2/var)\n",
        "print(\"组内相关：\",between_sd**2/var)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "B7059B12EE5B48EF854522F31D81D2FC",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 后验预测分布  \n",
        "\n",
        "* 可以看到相比于非池化模型，在层级模型中不同组的后验预测可信区间的长度、后验预测均值都更为接近"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 29,
      "metadata": {
        "collapsed": false,
        "id": "045A30C0EB3D441A91B29F132E284EBF",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Sampling: [y_est]\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "partial_ppc = pm.sample_posterior_predictive(partial_trace,\n",
        "                                            model=partial_pooled_model)\n",
        "partial_hdi_sum = ppc_sum(ppc=partial_ppc,\n",
        "                  data=df_first5)                                         "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 30,
      "metadata": {
        "collapsed": false,
        "id": "4E1317E216A248829E3742C3BB2A2710",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/040A9905AF2044FC9F0237E7B424C515/s5jqf5y45m.png\">"
            ],
            "text/plain": [
              "<Figure size 1500x600 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/040A9905AF2044FC9F0237E7B424C515/s5jqf9uh4a.png\">"
            ],
            "text/plain": [
              "<Figure size 1500x600 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/040A9905AF2044FC9F0237E7B424C515/s5jqfbvonj.png\">"
            ],
            "text/plain": [
              "<Figure size 1500x600 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "ppc_plot(hdi_sum=partial_hdi_sum)\n",
        "ppc_plot(hdi_sum=no_hdi_sum)\n",
        "ppc_plot(hdi_sum=complete_hdi_sum)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 31,
      "metadata": {
        "collapsed": false,
        "id": "55AA7A89A4634A4BB298256103819949",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>mean</th>\n",
              "      <th>sd</th>\n",
              "      <th>hdi_3%</th>\n",
              "      <th>hdi_97%</th>\n",
              "      <th>mcse_mean</th>\n",
              "      <th>mcse_sd</th>\n",
              "      <th>ess_bulk</th>\n",
              "      <th>ess_tail</th>\n",
              "      <th>r_hat</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>hyper_mu</th>\n",
              "      <td>40.629</td>\n",
              "      <td>0.729</td>\n",
              "      <td>39.337</td>\n",
              "      <td>42.100</td>\n",
              "      <td>0.007</td>\n",
              "      <td>0.005</td>\n",
              "      <td>10882.0</td>\n",
              "      <td>10954.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>mu[Kassel]</th>\n",
              "      <td>41.084</td>\n",
              "      <td>0.634</td>\n",
              "      <td>39.939</td>\n",
              "      <td>42.286</td>\n",
              "      <td>0.008</td>\n",
              "      <td>0.005</td>\n",
              "      <td>6816.0</td>\n",
              "      <td>8890.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>mu[METU]</th>\n",
              "      <td>40.055</td>\n",
              "      <td>0.513</td>\n",
              "      <td>39.070</td>\n",
              "      <td>41.009</td>\n",
              "      <td>0.005</td>\n",
              "      <td>0.004</td>\n",
              "      <td>9144.0</td>\n",
              "      <td>5299.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>mu[Oslo]</th>\n",
              "      <td>41.894</td>\n",
              "      <td>0.811</td>\n",
              "      <td>40.410</td>\n",
              "      <td>43.433</td>\n",
              "      <td>0.017</td>\n",
              "      <td>0.012</td>\n",
              "      <td>2231.0</td>\n",
              "      <td>742.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>mu[Southampton]</th>\n",
              "      <td>40.346</td>\n",
              "      <td>1.306</td>\n",
              "      <td>37.681</td>\n",
              "      <td>42.739</td>\n",
              "      <td>0.012</td>\n",
              "      <td>0.009</td>\n",
              "      <td>12110.0</td>\n",
              "      <td>10346.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>mu[Tsinghua]</th>\n",
              "      <td>39.773</td>\n",
              "      <td>0.544</td>\n",
              "      <td>38.722</td>\n",
              "      <td>40.757</td>\n",
              "      <td>0.005</td>\n",
              "      <td>0.004</td>\n",
              "      <td>10176.0</td>\n",
              "      <td>7878.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>within_variability</th>\n",
              "      <td>7.448</td>\n",
              "      <td>0.224</td>\n",
              "      <td>7.016</td>\n",
              "      <td>7.856</td>\n",
              "      <td>0.002</td>\n",
              "      <td>0.001</td>\n",
              "      <td>13778.0</td>\n",
              "      <td>11593.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>between_variability</th>\n",
              "      <td>1.211</td>\n",
              "      <td>0.635</td>\n",
              "      <td>0.237</td>\n",
              "      <td>2.344</td>\n",
              "      <td>0.012</td>\n",
              "      <td>0.008</td>\n",
              "      <td>1429.0</td>\n",
              "      <td>461.0</td>\n",
              "      <td>1.01</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "                       mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  \\\n",
              "hyper_mu             40.629  0.729  39.337   42.100      0.007    0.005   \n",
              "mu[Kassel]           41.084  0.634  39.939   42.286      0.008    0.005   \n",
              "mu[METU]             40.055  0.513  39.070   41.009      0.005    0.004   \n",
              "mu[Oslo]             41.894  0.811  40.410   43.433      0.017    0.012   \n",
              "mu[Southampton]      40.346  1.306  37.681   42.739      0.012    0.009   \n",
              "mu[Tsinghua]         39.773  0.544  38.722   40.757      0.005    0.004   \n",
              "within_variability    7.448  0.224   7.016    7.856      0.002    0.001   \n",
              "between_variability   1.211  0.635   0.237    2.344      0.012    0.008   \n",
              "\n",
              "                     ess_bulk  ess_tail  r_hat  \n",
              "hyper_mu              10882.0   10954.0   1.00  \n",
              "mu[Kassel]             6816.0    8890.0   1.00  \n",
              "mu[METU]               9144.0    5299.0   1.00  \n",
              "mu[Oslo]               2231.0     742.0   1.00  \n",
              "mu[Southampton]       12110.0   10346.0   1.00  \n",
              "mu[Tsinghua]          10176.0    7878.0   1.00  \n",
              "within_variability    13778.0   11593.0   1.00  \n",
              "between_variability    1429.0     461.0   1.01  "
            ]
          },
          "execution_count": 31,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "az.summary(partial_trace)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "0F809C7569F549CBA227AB8C99B4666C",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## Shrinkage & the bias-variance trade-off  \n",
        "\n",
        "从下图和表格中可以发现，三个模型的关系：  \n",
        "* 部分池化模型的超参数 (hyper_mu)的后验分布接近完全池化模型估计的参数 (mu)  \n",
        "* 部分池化模型中对站点参数的估计 (mu[?]) 的后验分布接近非池化模型估计的参数 (mu[?])  \n",
        "* 对于有的站点来说，相比于非池化模型，部分池化模型中的参数 (mu[?]) 更加靠近完全池化模型，这就是**分层模型的收缩 (shrinkage) 现象**。  \n",
        "\t* 例如，非池化模型中mu[Southampton]为38.6，而部分池化模型的mu[Southampton]=40.4，更加接近 完全池化模型的 mu 40.5。"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 32,
      "metadata": {
        "collapsed": false,
        "id": "61178FC197BE4613B52394CEDB92A615",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th>index</th>\n",
              "      <th>mean</th>\n",
              "      <th>sd</th>\n",
              "      <th>hdi_3%</th>\n",
              "      <th>hdi_97%</th>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>source</th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>Complete pool</th>\n",
              "      <th>0</th>\n",
              "      <td>mu</td>\n",
              "      <td>40.448</td>\n",
              "      <td>0.324</td>\n",
              "      <td>39.864</td>\n",
              "      <td>41.070</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th rowspan=\"5\" valign=\"top\">No pool</th>\n",
              "      <th>0</th>\n",
              "      <td>mu[Kassel]</td>\n",
              "      <td>41.342</td>\n",
              "      <td>0.732</td>\n",
              "      <td>39.927</td>\n",
              "      <td>42.662</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>mu[METU]</td>\n",
              "      <td>39.854</td>\n",
              "      <td>0.602</td>\n",
              "      <td>38.768</td>\n",
              "      <td>41.016</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>mu[Oslo]</td>\n",
              "      <td>42.699</td>\n",
              "      <td>0.748</td>\n",
              "      <td>41.249</td>\n",
              "      <td>44.059</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>mu[Southampton]</td>\n",
              "      <td>38.592</td>\n",
              "      <td>1.900</td>\n",
              "      <td>34.871</td>\n",
              "      <td>42.057</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>mu[Tsinghua]</td>\n",
              "      <td>39.474</td>\n",
              "      <td>0.449</td>\n",
              "      <td>38.639</td>\n",
              "      <td>40.288</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th rowspan=\"6\" valign=\"top\">Partial pool</th>\n",
              "      <th>0</th>\n",
              "      <td>hyper_mu</td>\n",
              "      <td>40.629</td>\n",
              "      <td>0.729</td>\n",
              "      <td>39.337</td>\n",
              "      <td>42.100</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>mu[Kassel]</td>\n",
              "      <td>41.084</td>\n",
              "      <td>0.634</td>\n",
              "      <td>39.939</td>\n",
              "      <td>42.286</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>mu[METU]</td>\n",
              "      <td>40.055</td>\n",
              "      <td>0.513</td>\n",
              "      <td>39.070</td>\n",
              "      <td>41.009</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>mu[Oslo]</td>\n",
              "      <td>41.894</td>\n",
              "      <td>0.811</td>\n",
              "      <td>40.410</td>\n",
              "      <td>43.433</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>mu[Southampton]</td>\n",
              "      <td>40.346</td>\n",
              "      <td>1.306</td>\n",
              "      <td>37.681</td>\n",
              "      <td>42.739</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>5</th>\n",
              "      <td>mu[Tsinghua]</td>\n",
              "      <td>39.773</td>\n",
              "      <td>0.544</td>\n",
              "      <td>38.722</td>\n",
              "      <td>40.757</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "                           index    mean     sd  hdi_3%  hdi_97%\n",
              "source                                                          \n",
              "Complete pool 0               mu  40.448  0.324  39.864   41.070\n",
              "No pool       0       mu[Kassel]  41.342  0.732  39.927   42.662\n",
              "              1         mu[METU]  39.854  0.602  38.768   41.016\n",
              "              2         mu[Oslo]  42.699  0.748  41.249   44.059\n",
              "              3  mu[Southampton]  38.592  1.900  34.871   42.057\n",
              "              4     mu[Tsinghua]  39.474  0.449  38.639   40.288\n",
              "Partial pool  0         hyper_mu  40.629  0.729  39.337   42.100\n",
              "              1       mu[Kassel]  41.084  0.634  39.939   42.286\n",
              "              2         mu[METU]  40.055  0.513  39.070   41.009\n",
              "              3         mu[Oslo]  41.894  0.811  40.410   43.433\n",
              "              4  mu[Southampton]  40.346  1.306  37.681   42.739\n",
              "              5     mu[Tsinghua]  39.773  0.544  38.722   40.757"
            ]
          },
          "execution_count": 32,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# 提取三个模型参数后验，筛选中其中含有mu参数的部分\n",
        "partial_stats = az.summary(partial_trace, var_names=[\"mu\"], filter_vars = \"like\", kind=\"stats\")\n",
        "no_stats = az.summary(no_trace, var_names=[\"mu\"], filter_vars = \"like\", kind=\"stats\")\n",
        "complete_stats = az.summary(complete_trace, var_names=[\"mu\"], filter_vars = \"like\", kind=\"stats\")\n",
        "# 设置一列，表明参数来源\n",
        "complete_stats['source'] = 'Complete pool'\n",
        "no_stats['source'] = 'No pool'\n",
        "partial_stats['source'] = 'Partial pool'\n",
        "# 合并三个模型的结果\n",
        "df_compare = pd.concat([complete_stats.reset_index(),\n",
        "                        no_stats.reset_index(),\n",
        "                        partial_stats.reset_index()])\n",
        "#设置索引，表明参数来源\n",
        "df_compare.set_index(['source', df_compare.index], inplace=True)\n",
        "df_compare"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 33,
      "metadata": {
        "collapsed": false,
        "id": "4A4E6350F5BE4C7E91D7799FF40BBA2B",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/F68ED1A3FD294BD797261CC378655645/s5jqfbll3d.png\">"
            ],
            "text/plain": [
              "<Figure size 2000x1000 with 3 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# 设置三个绘制坐标轴\n",
        "fig, (ax1, ax2, ax3) = plt.subplots(1,3, figsize=(20,10), sharex=True)\n",
        "\n",
        "# 绘制三个模型参数后验\n",
        "az.plot_forest(partial_trace, var_names=[\"mu\"], filter_vars = \"like\", combined=True, ax=ax1)\n",
        "ax1.set_title(\"Patial Pooling\")\n",
        "az.plot_forest(no_trace, var_names=[\"mu\"], filter_vars = \"like\", combined=True, ax=ax2)\n",
        "ax2.set_title(\"No Pooling\")\n",
        "az.plot_forest(complete_trace, var_names=[\"mu\"], filter_vars = \"like\", combined=True, ax=ax3)\n",
        "ax3.set_title(\"Complete Pooling\")\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "EB9A03E807BB43BD92463642A2C43FB2",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "* 在完全池化模型中，对于每一个站点来说，后验估计均值都是相同的，在弱先验的情况下，后验估计均值约等于所有观测值的平均值  \n",
        "$$  \n",
        "\\overline{y}_{\\text{global}} = \\frac{1}{n}\\sum_{\\text{all } i,j }y_{ij}  \n",
        "$$  \n",
        "\n",
        "* 在非池化模型中，对于每一个站点来说，其后验估计均值的估计只来自该组内部，在弱先验的情况下，后验估计均值约等于该组内所有观测值的平均值  \n",
        "$$  \n",
        "\\overline{y}_j = \\frac{1}{n_j}\\sum_{i=1}^{n_j} y_{ij}  \n",
        "$$  \n",
        "\n",
        "* 在层级模型中，后验估计均值则是在完全池化和非池化模型中找到一个平衡，组间(group-specific)参数有可能更倾向完全池化模型，也有可能更倾向于非池化模型，这种现象被称为shrinkage  \n",
        "> 当使用弱信息先验时，分层模型的后验平均预测结果（大致）是完全池化模型和非池化模型预测结果的加权平均  \n",
        "$$  \n",
        "\\frac{\\sigma^2_y}{\\sigma^2_y + n_j \\sigma^2_\\mu} \\overline{y}_{\\text{global}} + \\frac{n_j\\sigma^2_\\mu}{\\sigma^2_y + n_j \\sigma^2_\\mu} \\overline{y}_j  \n",
        "$$  \n",
        "\n",
        "* 我们可以分别绘制出每个观测值对应的后验预测均值，并观察这些后验预测均值在不同的模型之间发生了什么变化"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 34,
      "metadata": {
        "collapsed": false,
        "id": "6205851B39E847B496A59CB409C66126",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/1523843281AF468C951BA1836FCB390F/s5jqfczb1i.png\">"
            ],
            "text/plain": [
              "<Figure size 1500x500 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "#设置画布大小\n",
        "plt.figure(figsize=(15,5))\n",
        "#绘制完全池化模型下每个点对应的后验预测均值\n",
        "plt.scatter(complete_hdi_sum[\"obs_id\"],\n",
        "            complete_hdi_sum[\"mean\"],\n",
        "            alpha=0.15,\n",
        "            s=80,\n",
        "            label=\"Complete pooling\")\n",
        "#绘制非池化模型下每个点对应的后验预测均值\n",
        "plt.scatter(no_hdi_sum[\"obs_id\"],\n",
        "            no_hdi_sum[\"mean\"],\n",
        "            alpha=0.15,\n",
        "            s=80,\n",
        "            label=\"No pooling\")\n",
        "#绘制部分池化模型下每个点对应的后验预测均值\n",
        "plt.scatter(partial_hdi_sum[\"obs_id\"],\n",
        "            partial_hdi_sum[\"mean\"],\n",
        "            alpha=0.15,\n",
        "            s=80,\n",
        "            label=\"Partial pooling\")\n",
        "#设置y轴范围\n",
        "plt.ylim(38,44)\n",
        "#设置图例\n",
        "plt.legend()\n",
        "#计算每个站点的数据量，并根据数据量大小在x轴上进行刻度标识\n",
        "count_per_site = df_first5.groupby(level=\"Site\").size().values\n",
        "cumulative_count = count_per_site.cumsum()\n",
        "xtick = cumulative_count - count_per_site / 2\n",
        "plt.xticks(xtick,df_first5[\"Site\"].unique())\n",
        "#设置标题\n",
        "plt.title(\"Posterior mean of observed data\",\n",
        "          fontsize=15)\n",
        "sns.despine()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "7C7B48E399C84A9C9FD077B8E96A5F82",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "收缩是指在分层模型中，特定群体的局部趋势(各站点的估计结果)被拉向或收缩向全局趋势(所有数据的整体趋势)的现象。  \n",
        "\n",
        "$$  \n",
        "\\frac{\\sigma^2_y}{\\sigma^2_y + n_j \\sigma^2_\\mu} \\overline{y}_{\\text{global}} + \\frac{n_j\\sigma^2_\\mu}{\\sigma^2_y + n_j \\sigma^2_\\mu} \\overline{y}_j  \n",
        "$$  \n",
        "\n",
        "* 收缩的大小 (即完全池化模型均值和非池化模型均值的权重)，取决于站点的数量 $n_j$ 以及 组内和组间变异性（$\\sigma_y$和 $\\sigma_μ$）的比。  \n",
        "* 当站点的数量越少时，缩减率会增加。即我们越来越依赖全局趋势来了解一个组别很少 (可能意味着不可靠)的结果。  \n",
        "* 当组别内的变异性$\\sigma_y$与组别间的变异性 $\\sigma_μ$ 相比较大时，缩减会增加。即当各站点之间的差异很小时，我们会更依赖于用全局趋势来理解其中一个站点的数据。  \n",
        "> 这也说明，为什么 mu[Southampton] 受到的收缩影响最严重，因为 Southampton 站点的数据量 n 很少，并且组内 $\\sigma_y$ 变异很大。  \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "429131812F3F464F995F499D8E747D9E",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "层次模型在完全池化和无池化之间取得了平衡，优势在于：  \n",
        "\n",
        "1. 将对不同群体(站点)的观察结果推广到更广泛的总体中  \n",
        "2. 在了解任何单个群体时，借用所有其他群体的信息  \n",
        "\n",
        "这也导致了它的缺点：  \n",
        "* 如果站点数量少，并且组内变异大，那么我们对于总体和不同组别(站点)的估计都可能存在偏差"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "378E98DBBC3640AA90ECBB12BAFB53BD",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**bias-variance trade-off 偏差-方差权衡**  \n",
        "\n",
        "\n",
        "|               | 完全池化                       | 非池化                                 | 部分池化               |  \n",
        "| ------------- | ------------------------------ | -------------------------------------- | ---------------------- |  \n",
        "| 特点          | 仅考虑总体的变异，模型过于简单 | 仅考虑各组别的差异，结果难以推广到总体 | 同时考虑组间和组内变异 |  \n",
        "| 偏差-方差权衡 | 有较高的偏差和较低的方差       | 较低的偏差和较高的方差                 | 更为平衡               |  \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "1955519AE8E341939C597A8D2591B5DD",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## 预测新组的数据  \n",
        "\n",
        "* 我们可以根据当前的层级模型对新组别的数据进行预测，如\"Zurich\"站点  \n",
        "\n",
        "* 在pymc中，只要在`pm.sample_posterior_predictive`中传入层级模型的后验参数采样结果，即可以在层级模型的基础上对新数据生成预测  \n",
        "\n",
        "* 预测结果储存在`.predictions`中"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 35,
      "metadata": {
        "collapsed": false,
        "id": "9C1B547A2CA040BA9F50C557A0805BF4",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 选择站点为\"Zurich\"的数据\n",
        "new_group = df_raw[df_raw.Site==\"Zurich\"]\n",
        "# 生成被试索引\n",
        "new_group[\"obs_id\"] = range(len(new_group))\n",
        "# 生成站点索引\n",
        "new_group[\"site_idx\"] = pd.factorize(new_group.Site)[0]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 36,
      "metadata": {
        "collapsed": false,
        "id": "AF0B1A10E84746E29867995CE3BD4441",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Sampling: [new_mu, y_est]\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "new_coords = {\"site\": new_group[\"Site\"].unique(),\n",
        "          \"obs_id\": new_group.obs_id}\n",
        "\n",
        "with pm.Model(coords=new_coords) as partial_pooled_pred:\n",
        "    # Hyperpriors\n",
        "    var_y = pm.Exponential(\"within_variability\", 1)\n",
        "    var_mu = pm.Exponential(\"between_variability\", 1)\n",
        "    hyper_mu = pm.Normal(\"hyper_mu\", mu=40, sigma=20)\n",
        "    # 在这里我们需要传入一个新的参数名，因为传入的是一个新站点\n",
        "    new_mu = pm.Normal(\"new_mu\", mu=hyper_mu, sigma=var_mu, dims=\"site\")\n",
        "    # 其他的设置并没有发生改变\n",
        "    g = pm.MutableData(\"g\", new_group.site_idx, dims=\"obs_id\")\n",
        "\n",
        "    likelihood = pm.Normal(\"y_est\", mu=new_mu[g], sigma=var_y, observed=new_group.scontrol, dims=\"obs_id\")\n",
        "    # 进行后验预测估计，注意使用的是上一个模型的后验参数估计，partial_trace\n",
        "    pred_trace = pm.sample_posterior_predictive(partial_trace,\n",
        "                                                var_names=[\"y_est\"],\n",
        "                                                predictions=True,\n",
        "                                                extend_inferencedata=True,\n",
        "                                                random_seed=84735)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 37,
      "metadata": {
        "collapsed": false,
        "id": "CE5B7339E73547649EAA98D5B5653E68",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/0425E34030A94D30B7261166713A16C6/s5jqfdzbjn.png\">"
            ],
            "text/plain": [
              "<Figure size 1500x600 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "pred_hdi_sum = ppc_sum(ppc=pred_trace.predictions,data=new_group)\n",
        "ppc_plot(pred_hdi_sum)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "01DBF60E8F9B41E3AAA753B31CAB9886",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## bambi code  \n",
        "\n",
        "最后，我们演示如何通过 bambi 来实现三种模型。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "B9A7B3DB007148D7A18B93CA6FDCA17C",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**完全池化模型**  \n",
        "\n",
        "$$  \n",
        "\\begin{split}  \n",
        "Y_{ij} | \\mu, \\sigma & \\sim N(\\mu, \\sigma^2) \\\\  \n",
        "\\mu    & \\sim N(0, 50^2) \\\\  \n",
        "\\sigma & \\sim \\text{Exp}(1) \\\\  \n",
        "\\end{split}  \n",
        "$$  \n",
        "\n",
        "* 观测数据 Y 主要受到 $\\mu$ 的影响。在 bambi 中，可以通过 `\"scontrol ~ 1\"` 来表示这种关系。  \n",
        "* 结果中的 Intercept\t代表了  $\\mu$； scontrol_sigma 代表 $\\sigma$。"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 38,
      "metadata": {
        "collapsed": false,
        "id": "03144154A1ED4B7ABEFA69BE83E5814E",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Auto-assigning NUTS sampler...\n",
            "Initializing NUTS using jitter+adapt_diag...\n",
            "Multiprocess sampling (4 chains in 4 jobs)\n",
            "NUTS: [scontrol_sigma, Intercept]\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 2 seconds.\n"
          ]
        }
      ],
      "source": [
        "complete_model = bmb.Model(\"scontrol ~ 1\", df_first5,\n",
        "                           family=\"gaussian\")\n",
        "complete_idata = complete_model.fit(random_seed=84735)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 39,
      "metadata": {
        "collapsed": false,
        "id": "94D5BC2046D34CC3A0D58E8DE6AB6A6E",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/E5C145A97BC54221AA37E148A0B26F83/s5jqfmn6fp.svg\">"
            ],
            "text/plain": [
              "<graphviz.graphs.Digraph at 0x7f985bd7d7c0>"
            ]
          },
          "execution_count": 39,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "complete_model.graph()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 40,
      "metadata": {
        "collapsed": false,
        "id": "E234BCED326E49708D522D21F35C99A8",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>mean</th>\n",
              "      <th>sd</th>\n",
              "      <th>hdi_3%</th>\n",
              "      <th>hdi_97%</th>\n",
              "      <th>mcse_mean</th>\n",
              "      <th>mcse_sd</th>\n",
              "      <th>ess_bulk</th>\n",
              "      <th>ess_tail</th>\n",
              "      <th>r_hat</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>Intercept</th>\n",
              "      <td>40.452</td>\n",
              "      <td>0.328</td>\n",
              "      <td>39.854</td>\n",
              "      <td>41.064</td>\n",
              "      <td>0.005</td>\n",
              "      <td>0.004</td>\n",
              "      <td>3965.0</td>\n",
              "      <td>2954.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>scontrol_sigma</th>\n",
              "      <td>7.554</td>\n",
              "      <td>0.222</td>\n",
              "      <td>7.126</td>\n",
              "      <td>7.960</td>\n",
              "      <td>0.004</td>\n",
              "      <td>0.002</td>\n",
              "      <td>4002.0</td>\n",
              "      <td>2779.0</td>\n",
              "      <td>1.0</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "                  mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \\\n",
              "Intercept       40.452  0.328  39.854   41.064      0.005    0.004    3965.0   \n",
              "scontrol_sigma   7.554  0.222   7.126    7.960      0.004    0.002    4002.0   \n",
              "\n",
              "                ess_tail  r_hat  \n",
              "Intercept         2954.0    1.0  \n",
              "scontrol_sigma    2779.0    1.0  "
            ]
          },
          "execution_count": 40,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "az.summary(complete_idata)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "A1AFAE33385C4E9DA0520B74CF79B8D9",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**非池化模型**  \n",
        "\n",
        "$$  \n",
        "Y_{ij} | \\mu_j, \\sigma \\sim N(\\mu_j, \\sigma^2) \\\\  \n",
        "\n",
        "\\mu_j  \\sim N(0, 50^2) \\\\  \n",
        "\n",
        "\\sigma \\sim \\text{Exp}(1) \\\\  \n",
        "$$  \n",
        "\n",
        "* 该模型的关键在于让每一个站点拥有一个独特的参数 $\\mu_j$，对应 \"scontrol ~ 0 + C(Site)\"，其中的 0 表示去掉回归模型截距。  \n",
        "* 注意，不同于之前的 pymc 模型，bambi 默认认为 $\\sigma$ 在各站点间是相同的。  \n",
        "* 从结果可以看到：C(Site)[???] 代表不同站点的均值。"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 41,
      "metadata": {
        "collapsed": false,
        "id": "26B1B564EF424984AD99371118E41C58",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Auto-assigning NUTS sampler...\n",
            "Initializing NUTS using jitter+adapt_diag...\n",
            "Multiprocess sampling (4 chains in 4 jobs)\n",
            "NUTS: [scontrol_sigma, C(Site)]\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 4 seconds.\n"
          ]
        }
      ],
      "source": [
        "no_model = bmb.Model(\"scontrol ~ 0 + C(Site)\", df_first5,\n",
        "                     categorical=\"Site\")\n",
        "no_idata = no_model.fit(random_seed=84735)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 42,
      "metadata": {
        "collapsed": false,
        "id": "483EACA7049549E5A24955CA865AFE53",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/3CD4C40615B34907A36D507EF5A40BEC/s5jqg0h3ed.svg\">"
            ],
            "text/plain": [
              "<graphviz.graphs.Digraph at 0x7f985c83a100>"
            ]
          },
          "execution_count": 42,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "no_model.graph()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 43,
      "metadata": {
        "collapsed": false,
        "id": "AA4D3B7082094B3383CA0DC68E2E109D",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>mean</th>\n",
              "      <th>sd</th>\n",
              "      <th>hdi_3%</th>\n",
              "      <th>hdi_97%</th>\n",
              "      <th>mcse_mean</th>\n",
              "      <th>mcse_sd</th>\n",
              "      <th>ess_bulk</th>\n",
              "      <th>ess_tail</th>\n",
              "      <th>r_hat</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>C(Site)[Kassel]</th>\n",
              "      <td>41.297</td>\n",
              "      <td>0.749</td>\n",
              "      <td>39.910</td>\n",
              "      <td>42.723</td>\n",
              "      <td>0.010</td>\n",
              "      <td>0.007</td>\n",
              "      <td>5346.0</td>\n",
              "      <td>3268.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>C(Site)[METU]</th>\n",
              "      <td>39.856</td>\n",
              "      <td>0.543</td>\n",
              "      <td>38.804</td>\n",
              "      <td>40.868</td>\n",
              "      <td>0.008</td>\n",
              "      <td>0.005</td>\n",
              "      <td>5166.0</td>\n",
              "      <td>3211.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>C(Site)[Oslo]</th>\n",
              "      <td>42.684</td>\n",
              "      <td>0.808</td>\n",
              "      <td>41.197</td>\n",
              "      <td>44.210</td>\n",
              "      <td>0.010</td>\n",
              "      <td>0.007</td>\n",
              "      <td>6134.0</td>\n",
              "      <td>3793.0</td>\n",
              "      <td>1.01</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>C(Site)[Southampton]</th>\n",
              "      <td>38.677</td>\n",
              "      <td>3.094</td>\n",
              "      <td>33.393</td>\n",
              "      <td>45.077</td>\n",
              "      <td>0.041</td>\n",
              "      <td>0.029</td>\n",
              "      <td>5751.0</td>\n",
              "      <td>3327.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>C(Site)[Tsinghua]</th>\n",
              "      <td>39.468</td>\n",
              "      <td>0.567</td>\n",
              "      <td>38.407</td>\n",
              "      <td>40.492</td>\n",
              "      <td>0.007</td>\n",
              "      <td>0.005</td>\n",
              "      <td>6396.0</td>\n",
              "      <td>3458.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>scontrol_sigma</th>\n",
              "      <td>7.491</td>\n",
              "      <td>0.227</td>\n",
              "      <td>7.065</td>\n",
              "      <td>7.912</td>\n",
              "      <td>0.003</td>\n",
              "      <td>0.002</td>\n",
              "      <td>5506.0</td>\n",
              "      <td>3170.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "                        mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  \\\n",
              "C(Site)[Kassel]       41.297  0.749  39.910   42.723      0.010    0.007   \n",
              "C(Site)[METU]         39.856  0.543  38.804   40.868      0.008    0.005   \n",
              "C(Site)[Oslo]         42.684  0.808  41.197   44.210      0.010    0.007   \n",
              "C(Site)[Southampton]  38.677  3.094  33.393   45.077      0.041    0.029   \n",
              "C(Site)[Tsinghua]     39.468  0.567  38.407   40.492      0.007    0.005   \n",
              "scontrol_sigma         7.491  0.227   7.065    7.912      0.003    0.002   \n",
              "\n",
              "                      ess_bulk  ess_tail  r_hat  \n",
              "C(Site)[Kassel]         5346.0    3268.0   1.00  \n",
              "C(Site)[METU]           5166.0    3211.0   1.00  \n",
              "C(Site)[Oslo]           6134.0    3793.0   1.01  \n",
              "C(Site)[Southampton]    5751.0    3327.0   1.00  \n",
              "C(Site)[Tsinghua]       6396.0    3458.0   1.00  \n",
              "scontrol_sigma          5506.0    3170.0   1.00  "
            ]
          },
          "execution_count": 43,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "az.summary(no_idata)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "6C130F06F243475DAEF23F78BC08272D",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**部分池化模型**  \n",
        "\n",
        "$$  \n",
        "\\begin{split}  \n",
        "Y_{ij} | \\mu_j, \\sigma_y & \\sim N(\\mu_j, \\sigma_y^2) \\;\\; \\text{ with } \\;\\; \\mu_j = \\mu + b_{j}  \\\\  \n",
        "b_{j} | \\sigma_\\mu    & \\stackrel{ind}{\\sim} N(0, \\sigma_\\mu^2) \\\\  \n",
        "\\mu           & \\sim N(40, 20^2) \\\\  \n",
        "\\sigma_y      & \\sim \\text{Exp}(1) \\\\  \n",
        "\\sigma_\\mu    & \\sim \\text{Exp}(1) \\\\  \n",
        "\\end{split}  \n",
        "$$  \n",
        "\n",
        "* 部分池化模型的关键在于，参数的层级依赖关系。  \n",
        "* 而 bambi 提供了更简单的方法构建分层模型 `\"scontrol ~ 1 + (1|Site)\"`。  \n",
        "\t*  在完全池化模型 `\"scontrol ~ 1\"` 的基础上加入 `(1|Site)` 分层结构，就定义了好了分层模型。  \n",
        "\t*  注意，bambi 默认使用之前提到的第二种分层模型的定义形式 (线性模型的定义形式)  \n",
        "*  从结果可以看到：  \n",
        "\t*  Intercept 代表了全局总体的均值，1|Site_sigma 代表了 不同站点形成的变异，即组间变异。  \n",
        "\t*  1|Site[???] 代表不同站点相对于全局总体的均值 Intercept 的偏移，例如，站点 Kassel\t的均值 = Intercept + 1|Site[Kassel]\t= 40.61+0.57 = 41.18  \n",
        "\t*  scontrol_sigma 为组内变异。"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 44,
      "metadata": {
        "collapsed": false,
        "id": "73A11CC0E91D4C678E925908AB0E6064",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Auto-assigning NUTS sampler...\n",
            "Initializing NUTS using jitter+adapt_diag...\n",
            "Multiprocess sampling (4 chains in 4 jobs)\n",
            "NUTS: [scontrol_sigma, Intercept, 1|Site_sigma, 1|Site_offset]\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "\n",
              "<style>\n",
              "    /* Turns off some styling */\n",
              "    progress {\n",
              "        /* gets rid of default border in Firefox and Opera. */\n",
              "        border: none;\n",
              "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
              "        background-size: auto;\n",
              "    }\n",
              "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
              "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
              "    }\n",
              "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
              "        background: #F44336;\n",
              "    }\n",
              "</style>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 9 seconds.\n",
            "There were 6 divergences after tuning. Increase `target_accept` or reparameterize.\n"
          ]
        }
      ],
      "source": [
        "partial_model = bmb.Model(\"scontrol ~ 1 + (1|Site)\", df_first5,\n",
        "                     categorical=\"Site\")\n",
        "partial_idata = partial_model.fit(random_seed=84735)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 45,
      "metadata": {
        "collapsed": false,
        "id": "1C091EF585964BFAA433ADE8AE1E312A",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<img src=\"https://cdn.kesci.com/upload/rt/EA25CE759B8A47418252A8C2E40C601F/s5jqgoamag.svg\">"
            ],
            "text/plain": [
              "<graphviz.graphs.Digraph at 0x7f985bfea340>"
            ]
          },
          "execution_count": 45,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "partial_model.graph()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 46,
      "metadata": {
        "collapsed": false,
        "id": "3B46789CCED9494394F3A8B4EB12DBF0",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "scrolled": false,
        "slideshow": {
          "slide_type": "fragment"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>mean</th>\n",
              "      <th>sd</th>\n",
              "      <th>hdi_3%</th>\n",
              "      <th>hdi_97%</th>\n",
              "      <th>mcse_mean</th>\n",
              "      <th>mcse_sd</th>\n",
              "      <th>ess_bulk</th>\n",
              "      <th>ess_tail</th>\n",
              "      <th>r_hat</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>Intercept</th>\n",
              "      <td>40.612</td>\n",
              "      <td>1.009</td>\n",
              "      <td>38.780</td>\n",
              "      <td>42.748</td>\n",
              "      <td>0.054</td>\n",
              "      <td>0.039</td>\n",
              "      <td>439.0</td>\n",
              "      <td>376.0</td>\n",
              "      <td>1.01</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>scontrol_sigma</th>\n",
              "      <td>7.489</td>\n",
              "      <td>0.225</td>\n",
              "      <td>7.084</td>\n",
              "      <td>7.912</td>\n",
              "      <td>0.004</td>\n",
              "      <td>0.003</td>\n",
              "      <td>2907.0</td>\n",
              "      <td>2551.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1|Site_sigma</th>\n",
              "      <td>1.909</td>\n",
              "      <td>1.123</td>\n",
              "      <td>0.134</td>\n",
              "      <td>4.068</td>\n",
              "      <td>0.043</td>\n",
              "      <td>0.031</td>\n",
              "      <td>750.0</td>\n",
              "      <td>705.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1|Site[Kassel]</th>\n",
              "      <td>0.575</td>\n",
              "      <td>1.112</td>\n",
              "      <td>-1.665</td>\n",
              "      <td>2.750</td>\n",
              "      <td>0.053</td>\n",
              "      <td>0.047</td>\n",
              "      <td>532.0</td>\n",
              "      <td>394.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1|Site[METU]</th>\n",
              "      <td>-0.629</td>\n",
              "      <td>1.065</td>\n",
              "      <td>-2.897</td>\n",
              "      <td>1.274</td>\n",
              "      <td>0.054</td>\n",
              "      <td>0.039</td>\n",
              "      <td>453.0</td>\n",
              "      <td>381.0</td>\n",
              "      <td>1.01</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1|Site[Oslo]</th>\n",
              "      <td>1.573</td>\n",
              "      <td>1.195</td>\n",
              "      <td>-0.537</td>\n",
              "      <td>3.846</td>\n",
              "      <td>0.053</td>\n",
              "      <td>0.044</td>\n",
              "      <td>603.0</td>\n",
              "      <td>453.0</td>\n",
              "      <td>1.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1|Site[Southampton]</th>\n",
              "      <td>-0.467</td>\n",
              "      <td>1.748</td>\n",
              "      <td>-3.973</td>\n",
              "      <td>2.636</td>\n",
              "      <td>0.070</td>\n",
              "      <td>0.049</td>\n",
              "      <td>807.0</td>\n",
              "      <td>592.0</td>\n",
              "      <td>1.01</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1|Site[Tsinghua]</th>\n",
              "      <td>-0.946</td>\n",
              "      <td>1.089</td>\n",
              "      <td>-3.291</td>\n",
              "      <td>0.925</td>\n",
              "      <td>0.052</td>\n",
              "      <td>0.037</td>\n",
              "      <td>518.0</td>\n",
              "      <td>411.0</td>\n",
              "      <td>1.01</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "                       mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  \\\n",
              "Intercept            40.612  1.009  38.780   42.748      0.054    0.039   \n",
              "scontrol_sigma        7.489  0.225   7.084    7.912      0.004    0.003   \n",
              "1|Site_sigma          1.909  1.123   0.134    4.068      0.043    0.031   \n",
              "1|Site[Kassel]        0.575  1.112  -1.665    2.750      0.053    0.047   \n",
              "1|Site[METU]         -0.629  1.065  -2.897    1.274      0.054    0.039   \n",
              "1|Site[Oslo]          1.573  1.195  -0.537    3.846      0.053    0.044   \n",
              "1|Site[Southampton]  -0.467  1.748  -3.973    2.636      0.070    0.049   \n",
              "1|Site[Tsinghua]     -0.946  1.089  -3.291    0.925      0.052    0.037   \n",
              "\n",
              "                     ess_bulk  ess_tail  r_hat  \n",
              "Intercept               439.0     376.0   1.01  \n",
              "scontrol_sigma         2907.0    2551.0   1.00  \n",
              "1|Site_sigma            750.0     705.0   1.00  \n",
              "1|Site[Kassel]          532.0     394.0   1.00  \n",
              "1|Site[METU]            453.0     381.0   1.01  \n",
              "1|Site[Oslo]            603.0     453.0   1.00  \n",
              "1|Site[Southampton]     807.0     592.0   1.01  \n",
              "1|Site[Tsinghua]        518.0     411.0   1.01  "
            ]
          },
          "execution_count": 46,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "az.summary(partial_idata)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "E77877F9FF6348FABD585A77445F37D9",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**🎯开放式练习**  \n",
        "\n",
        "尝试使用 bambi 构建分层模型：  \n",
        "1. 自行选取其他站点的数据进行分析。  \n",
        "2. 可以尝试使用部分池化中的第一种公式构建分层模型，提示：可使用`\"scontrol ~ 0 + (1|Site)\"`构建模型。  \n",
        "3. 比较组间组内变异的差异。  \n",
        "4. 使用 `az.plot_forest` 绘制站点参数森林图。  \n",
        "5. 思考分层模型带来的收缩效应 (shrinkake)。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "E872B08076A94972B41792422D579F73",
        "jupyter": {},
        "notebookId": "6578310fc5a8cddbdabda6f8",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## 总结  \n",
        "\n",
        "在本节课中，我们学习了部分池化和层级模型(hierarchical model)：  \n",
        "- 了解层级数据结构的形式，区分了三种建模方法：完全池化(complete pooling),  非池化(no pooling), 部分池化(patial pooling)  \n",
        "- 讨论了组间变异(between variability)和组内变异(within variability)的差异，以及偏差方差权衡  \n",
        "- 通过pymc 和 bambi 实现三种不同的模型  \n"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.11.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}
